157 lines
5.2 KiB
C
157 lines
5.2 KiB
C
/*
|
|
* Single-precision vector e^x function.
|
|
*
|
|
* Copyright (c) 2019-2023, Arm Limited.
|
|
* SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
|
|
*/
|
|
|
|
#include "sv_math.h"
|
|
#include "pl_sig.h"
|
|
#include "pl_test.h"
|
|
|
|
#if SV_SUPPORTED
|
|
|
|
#define C(i) __sv_expf_poly[i]
|
|
|
|
#define InvLn2 (0x1.715476p+0f)
|
|
#define Ln2hi (0x1.62e4p-1f)
|
|
#define Ln2lo (0x1.7f7d1cp-20f)
|
|
|
|
#if SV_EXPF_USE_FEXPA
|
|
|
|
#define Shift (0x1.903f8p17f) /* 1.5*2^17 + 127. */
|
|
#define Thres \
|
|
(0x1.5d5e2ap+6f) /* Roughly 87.3. For x < -Thres, the result is subnormal \
|
|
and not handled correctly by FEXPA. */
|
|
|
|
static NOINLINE sv_f32_t
|
|
special_case (sv_f32_t x, sv_f32_t y, svbool_t special)
|
|
{
|
|
/* The special-case handler from the Neon routine does not handle subnormals
|
|
in a way that is compatible with FEXPA. For the FEXPA variant we just fall
|
|
back to scalar expf. */
|
|
return sv_call_f32 (expf, x, y, special);
|
|
}
|
|
|
|
#else
|
|
|
|
#define Shift (0x1.8p23f) /* 1.5 * 2^23. */
|
|
#define Thres (126.0f)
|
|
|
|
/* Special-case handler adapted from Neon variant. Uses s, y and n to produce
|
|
the final result (normal cases included). It performs an update of all lanes!
|
|
Therefore:
|
|
- all previous computation need to be done on all lanes indicated by input
|
|
pg
|
|
- we cannot simply apply the special case to the special-case-activated
|
|
lanes. Besides it is likely that this would not increase performance (no
|
|
scatter/gather). */
|
|
static inline sv_f32_t
|
|
specialcase (svbool_t pg, sv_f32_t poly, sv_f32_t n, sv_u32_t e,
|
|
svbool_t p_cmp1, sv_f32_t scale)
|
|
{
|
|
/* s=2^(n/N) may overflow, break it up into s=s1*s2,
|
|
such that exp = s + s*y can be computed as s1*(s2+s2*y)
|
|
and s1*s1 overflows only if n>0. */
|
|
|
|
/* If n<=0 then set b to 0x820...0, 0 otherwise. */
|
|
svbool_t p_sign = svcmple_n_f32 (pg, n, 0.0f); /* n <= 0. */
|
|
sv_u32_t b
|
|
= svdup_n_u32_z (p_sign, 0x82000000); /* Inactive lanes set to 0. */
|
|
|
|
/* Set s1 to generate overflow depending on sign of exponent n. */
|
|
sv_f32_t s1
|
|
= sv_as_f32_u32 (svadd_n_u32_x (pg, b, 0x7f000000)); /* b + 0x7f000000. */
|
|
/* Offset s to avoid overflow in final result if n is below threshold. */
|
|
sv_f32_t s2 = sv_as_f32_u32 (
|
|
svsub_u32_x (pg, e, b)); /* as_u32 (s) - 0x3010...0 + b. */
|
|
|
|
/* |n| > 192 => 2^(n/N) overflows. */
|
|
svbool_t p_cmp2 = svacgt_n_f32 (pg, n, 192.0f);
|
|
|
|
sv_f32_t r2 = svmul_f32_x (pg, s1, s1);
|
|
sv_f32_t r1 = sv_fma_f32_x (pg, poly, s2, s2);
|
|
r1 = svmul_f32_x (pg, r1, s1);
|
|
sv_f32_t r0 = sv_fma_f32_x (pg, poly, scale, scale);
|
|
|
|
/* Apply condition 1 then 2.
|
|
Returns r2 if cond2 is true, otherwise
|
|
if cond1 is true then return r1, otherwise return r0. */
|
|
sv_f32_t r = svsel_f32 (p_cmp1, r1, r0);
|
|
|
|
return svsel_f32 (p_cmp2, r2, r);
|
|
}
|
|
|
|
#endif
|
|
|
|
/* Optimised single-precision SVE exp function. By default this is an SVE port
|
|
of the Neon algorithm from math/. Alternatively, enable a modification of
|
|
that algorithm that looks up scale using SVE FEXPA instruction with
|
|
SV_EXPF_USE_FEXPA.
|
|
|
|
Worst-case error of the default algorithm is 1.95 ulp:
|
|
__sv_expf(-0x1.4cb74ap+2) got 0x1.6a022cp-8
|
|
want 0x1.6a023p-8.
|
|
|
|
Worst-case error when using FEXPA is 1.04 ulp:
|
|
__sv_expf(0x1.a8eda4p+1) got 0x1.ba74bcp+4
|
|
want 0x1.ba74bap+4. */
|
|
sv_f32_t
|
|
__sv_expf_x (sv_f32_t x, const svbool_t pg)
|
|
{
|
|
/* exp(x) = 2^n (1 + poly(r)), with 1 + poly(r) in [1/sqrt(2),sqrt(2)]
|
|
x = ln2*n + r, with r in [-ln2/2, ln2/2]. */
|
|
|
|
/* n = round(x/(ln2/N)). */
|
|
sv_f32_t z = sv_fma_n_f32_x (pg, InvLn2, x, sv_f32 (Shift));
|
|
sv_f32_t n = svsub_n_f32_x (pg, z, Shift);
|
|
|
|
/* r = x - n*ln2/N. */
|
|
sv_f32_t r = sv_fma_n_f32_x (pg, -Ln2hi, n, x);
|
|
r = sv_fma_n_f32_x (pg, -Ln2lo, n, r);
|
|
|
|
/* scale = 2^(n/N). */
|
|
#if SV_EXPF_USE_FEXPA
|
|
/* NaNs also need special handling with FEXPA. */
|
|
svbool_t is_special_case
|
|
= svorr_b_z (pg, svacgt_n_f32 (pg, x, Thres), svcmpne_f32 (pg, x, x));
|
|
sv_f32_t scale = svexpa_f32 (sv_as_u32_f32 (z));
|
|
#else
|
|
sv_u32_t e = svlsl_n_u32_x (pg, sv_as_u32_f32 (z), 23);
|
|
svbool_t is_special_case = svacgt_n_f32 (pg, n, Thres);
|
|
sv_f32_t scale = sv_as_f32_u32 (svadd_n_u32_x (pg, e, 0x3f800000));
|
|
#endif
|
|
|
|
/* y = exp(r) - 1 ~= r + C1 r^2 + C2 r^3 + C3 r^4. */
|
|
sv_f32_t r2 = svmul_f32_x (pg, r, r);
|
|
sv_f32_t p = sv_fma_n_f32_x (pg, C (0), r, sv_f32 (C (1)));
|
|
sv_f32_t q = sv_fma_n_f32_x (pg, C (2), r, sv_f32 (C (3)));
|
|
q = sv_fma_f32_x (pg, p, r2, q);
|
|
p = svmul_n_f32_x (pg, r, C (4));
|
|
sv_f32_t poly = sv_fma_f32_x (pg, q, r2, p);
|
|
|
|
if (unlikely (svptest_any (pg, is_special_case)))
|
|
#if SV_EXPF_USE_FEXPA
|
|
return special_case (x, sv_fma_f32_x (pg, poly, scale, scale),
|
|
is_special_case);
|
|
#else
|
|
return specialcase (pg, poly, n, e, is_special_case, scale);
|
|
#endif
|
|
|
|
return sv_fma_f32_x (pg, poly, scale, scale);
|
|
}
|
|
|
|
PL_ALIAS (__sv_expf_x, _ZGVsMxv_expf)
|
|
|
|
PL_SIG (SV, F, 1, exp, -9.9, 9.9)
|
|
PL_TEST_ULP (__sv_expf, 1.46)
|
|
PL_TEST_INTERVAL (__sv_expf, 0, 0x1p-23, 40000)
|
|
PL_TEST_INTERVAL (__sv_expf, 0x1p-23, 1, 50000)
|
|
PL_TEST_INTERVAL (__sv_expf, 1, 0x1p23, 50000)
|
|
PL_TEST_INTERVAL (__sv_expf, 0x1p23, inf, 50000)
|
|
PL_TEST_INTERVAL (__sv_expf, -0, -0x1p-23, 40000)
|
|
PL_TEST_INTERVAL (__sv_expf, -0x1p-23, -1, 50000)
|
|
PL_TEST_INTERVAL (__sv_expf, -1, -0x1p23, 50000)
|
|
PL_TEST_INTERVAL (__sv_expf, -0x1p23, -inf, 50000)
|
|
#endif // SV_SUPPORTED
|