diff options
author | Joe Ramsay <Joe.Ramsay@arm.com> | 2023-01-19 13:19:38 +0000 |
---|---|---|
committer | Joe Ramsay <joe.ramsay@arm.com> | 2023-01-19 13:19:43 +0000 |
commit | 76c2badba375dd8f785e92079dbd4290038a2750 (patch) | |
tree | 64325bba339d4ac084d0b3b46a3e85e99b1400ad | |
parent | cd28f3c4253257027f9834eb7b17db047346a267 (diff) | |
download | arm-optimized-routines-76c2badba375dd8f785e92079dbd4290038a2750.tar.gz |
pl/math: Add vector/Neon tan
New routine uses a similar technique to the single-precision Neon
routine, but with an extra reduction to pi/8 using the double-angle
formula. It is accurate to 3.5 ULP.
-rw-r--r-- | pl/math/include/mathlib.h | 4 | ||||
-rw-r--r-- | pl/math/math_config.h | 5 | ||||
-rw-r--r-- | pl/math/s_tan_3u5.c | 6 | ||||
-rw-r--r-- | pl/math/tools/tan.sollya | 20 | ||||
-rw-r--r-- | pl/math/v_tan_3u5.c | 102 | ||||
-rw-r--r-- | pl/math/v_tan_data.c | 15 | ||||
-rw-r--r-- | pl/math/vn_tan_3u5.c | 12 |
7 files changed, 164 insertions, 0 deletions
diff --git a/pl/math/include/mathlib.h b/pl/math/include/mathlib.h index 74c3b34..af5f9f9 100644 --- a/pl/math/include/mathlib.h +++ b/pl/math/include/mathlib.h @@ -70,6 +70,7 @@ double __s_log10 (double); double __s_log1p (double); double __s_log2 (double); double __s_sinh (double); +double __s_tan (double); double __s_tanh (double); #if __aarch64__ @@ -113,6 +114,7 @@ __f64x2_t __v_log2 (__f64x2_t); __f32x4_t __v_sinhf (__f32x4_t); __f64x2_t __v_sinh (__f64x2_t); __f32x4_t __v_tanf (__f32x4_t); +__f64x2_t __v_tan (__f64x2_t); __f32x4_t __v_tanhf (__f32x4_t); __f64x2_t __v_tanh (__f64x2_t); @@ -149,6 +151,7 @@ __vpcs __f64x2_t __vn_log2 (__f64x2_t); __vpcs __f32x4_t __vn_sinhf (__f32x4_t); __vpcs __f64x2_t __vn_sinh (__f64x2_t); __vpcs __f32x4_t __vn_tanf (__f32x4_t); +__vpcs __f64x2_t __vn_tan (__f64x2_t); __vpcs __f32x4_t __vn_tanhf (__f32x4_t); __vpcs __f64x2_t __vn_tanh (__f64x2_t); @@ -182,6 +185,7 @@ __vpcs __f64x2_t _ZGVnN2v_log2 (__f64x2_t); __vpcs __f32x4_t _ZGVnN4v_sinhf (__f32x4_t); __vpcs __f64x2_t _ZGVnN2v_sinh (__f64x2_t); __vpcs __f32x4_t _ZGVnN4v_tanf (__f32x4_t); +__vpcs __f64x2_t _ZGVnN2v_tan (__f64x2_t); __vpcs __f32x4_t _ZGVnN4v_tanhf (__f32x4_t); __vpcs __f64x2_t _ZGVnN2v_tanh (__f64x2_t); diff --git a/pl/math/math_config.h b/pl/math/math_config.h index 03b4ad4..9a7ce96 100644 --- a/pl/math/math_config.h +++ b/pl/math/math_config.h @@ -564,4 +564,9 @@ extern const struct cbrt_data double table[5]; } __cbrt_data HIDDEN; +extern const struct v_tan_data +{ + double neg_half_pi_hi, neg_half_pi_lo; + double poly[9]; +} __v_tan_data HIDDEN; #endif diff --git a/pl/math/s_tan_3u5.c b/pl/math/s_tan_3u5.c new file mode 100644 index 0000000..adb807c --- /dev/null +++ b/pl/math/s_tan_3u5.c @@ -0,0 +1,6 @@ +/* + * Copyright (c) 2023, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ +#define SCALAR 1 +#include "v_tan_3u5.c" diff --git a/pl/math/tools/tan.sollya b/pl/math/tools/tan.sollya new file mode 100644 index 0000000..bb0bb28 --- /dev/null +++ b/pl/math/tools/tan.sollya @@ -0,0 +1,20 @@ +// polynomial for approximating double precision tan(x) +// +// Copyright (c) 2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 8; + +// interval bounds +a = 0x1.0p-126; +b = pi / 8; + +display = hexadecimal; + +f = (tan(sqrt(x))-sqrt(x))/x^(3/2); +poly = fpminimax(f, deg, [|double ...|], [a*a;b*b]); + +//print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30)); +print("in [",a,b,"]"); +print("coeffs:"); +for i from 0 to deg do coeff(poly,i); diff --git a/pl/math/v_tan_3u5.c b/pl/math/v_tan_3u5.c new file mode 100644 index 0000000..f87bacc --- /dev/null +++ b/pl/math/v_tan_3u5.c @@ -0,0 +1,102 @@ +/* + * Double-precision vector tan(x) function. + * + * Copyright (c) 2023, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ + +#include "v_math.h" +#include "estrin.h" +#include "pl_sig.h" +#include "pl_test.h" + +#if V_SUPPORTED + +#define MHalfPiHi v_f64 (__v_tan_data.neg_half_pi_hi) +#define MHalfPiLo v_f64 (__v_tan_data.neg_half_pi_lo) +#define TwoOverPi v_f64 (0x1.45f306dc9c883p-1) +#define Shift v_f64 (0x1.8p52) +#define AbsMask 0x7fffffffffffffff +#define RangeVal 0x4160000000000000 /* asuint64(2^23). */ +#define TinyBound 0x3e50000000000000 /* asuint64(2^-26). */ +#define C(i) v_f64 (__v_tan_data.poly[i]) + +/* Special cases (fall back to scalar calls). */ +VPCS_ATTR +NOINLINE static v_f64_t +specialcase (v_f64_t x) +{ + return v_call_f64 (tan, x, x, v_u64 (-1)); +} + +/* Vector approximation for double-precision tan. + Maximum measured error is 3.48 ULP: + __v_tan(0x1.4457047ef78d8p+20) got -0x1.f6ccd8ecf7dedp+37 + want -0x1.f6ccd8ecf7deap+37. */ +VPCS_ATTR +v_f64_t V_NAME (tan) (v_f64_t x) +{ + v_u64_t iax = v_as_u64_f64 (x) & AbsMask; + + /* Our argument reduction cannot calculate q with sufficient accuracy for very + large inputs. Fall back to scalar routine for all lanes if any are too + large, or Inf/NaN. If fenv exceptions are expected, also fall back for tiny + input to avoid underflow. Note pl does not supply a scalar double-precision + tan, so the fallback will be statically linked from the system libm. */ +#if WANT_SIMD_EXCEPT + if (unlikely (v_any_u64 (iax - TinyBound > RangeVal - TinyBound))) +#else + if (unlikely (v_any_u64 (iax > RangeVal))) +#endif + return specialcase (x); + + /* q = nearest integer to 2 * x / pi. */ + v_f64_t q = v_fma_f64 (x, TwoOverPi, Shift) - Shift; + v_s64_t qi = v_to_s64_f64 (q); + + /* Use q to reduce x to r in [-pi/4, pi/4], by: + r = x - q * pi/2, in extended precision. */ + v_f64_t r = x; + r = v_fma_f64 (q, MHalfPiHi, r); + r = v_fma_f64 (q, MHalfPiLo, r); + /* Further reduce r to [-pi/8, pi/8], to be reconstructed using double angle + formula. */ + r = r * 0.5; + + /* Approximate tan(r) using order 8 polynomial. + tan(x) is odd, so polynomial has the form: + tan(x) ~= x + C0 * x^3 + C1 * x^5 + C3 * x^7 + ... + Hence we first approximate P(r) = C1 + C2 * r^2 + C3 * r^4 + ... + Then compute the approximation by: + tan(r) ~= r + r^3 * (C0 + r^2 * P(r)). */ + v_f64_t r2 = r * r, r4 = r2 * r2, r8 = r4 * r4; + /* Use offset version of Estrin wrapper to evaluate from C1 onwards. */ + v_f64_t p = ESTRIN_7_ (r2, r4, r8, C, 1); + p = v_fma_f64 (p, r2, C (0)); + p = v_fma_f64 (r2, p * r, r); + + /* Recombination uses double-angle formula: + tan(2x) = 2 * tan(x) / (1 - (tan(x))^2) + and reciprocity around pi/2: + tan(x) = 1 / (tan(pi/2 - x)) + to assemble result using change-of-sign and conditional selection of + numerator/denominator, dependent on odd/even-ness of q (hence quadrant). */ + v_f64_t n = v_fma_f64 (p, p, v_f64 (-1)); + v_f64_t d = p * 2; + + v_u64_t use_recip = v_cond_u64 ((v_as_u64_s64 (qi) & 1) == 0); + + return v_sel_f64 (use_recip, -d, n) / v_sel_f64 (use_recip, n, d); +} +VPCS_ALIAS + +PL_SIG (V, D, 1, tan, -3.1, 3.1) +PL_TEST_ULP (V_NAME (tan), 2.99) +PL_TEST_EXPECT_FENV (V_NAME (tan), WANT_SIMD_EXCEPT) +PL_TEST_INTERVAL (V_NAME (tan), 0, TinyBound, 5000) +PL_TEST_INTERVAL (V_NAME (tan), TinyBound, RangeVal, 100000) +PL_TEST_INTERVAL (V_NAME (tan), RangeVal, inf, 5000) +PL_TEST_INTERVAL (V_NAME (tan), -0, -TinyBound, 5000) +PL_TEST_INTERVAL (V_NAME (tan), -TinyBound, -RangeVal, 100000) +PL_TEST_INTERVAL (V_NAME (tan), -RangeVal, -inf, 5000) +#endif diff --git a/pl/math/v_tan_data.c b/pl/math/v_tan_data.c new file mode 100644 index 0000000..04e2516 --- /dev/null +++ b/pl/math/v_tan_data.c @@ -0,0 +1,15 @@ +/* + * Coefficients and helpers for double-precision vector tan(x) function. + * + * Copyright (c) 2023, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ +#include "math_config.h" + +const struct v_tan_data __v_tan_data + = {.neg_half_pi_hi = -0x1.921fb54442d18p0, + .neg_half_pi_lo = -0x1.1a62633145c07p-54, + .poly + = {0x1.5555555555556p-2, 0x1.1111111110a63p-3, 0x1.ba1ba1bb46414p-5, + 0x1.664f47e5b5445p-6, 0x1.226e5e5ecdfa3p-7, 0x1.d6c7ddbf87047p-9, + 0x1.7ea75d05b583ep-10, 0x1.289f22964a03cp-11, 0x1.4e4fd14147622p-12}}; diff --git a/pl/math/vn_tan_3u5.c b/pl/math/vn_tan_3u5.c new file mode 100644 index 0000000..a4efb06 --- /dev/null +++ b/pl/math/vn_tan_3u5.c @@ -0,0 +1,12 @@ +/* + * AdvSIMD vector PCS variant of __v_tan. + * + * Copyright (c) 2023, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ +#include "include/mathlib.h" +#ifdef __vpcs +#define VPCS 1 +#define VPCS_ALIAS PL_ALIAS (__vn_tan, _ZGVnN2v_tan) +#include "v_tan_3u5.c" +#endif |