mirror of
git://sourceware.org/git/glibc.git
synced 2026-09-08 23:58:31 +08:00
Pierre Blanchard
2ff2565df7
Vector variants of the new C23 powr routines. These provide same maximum error error as pow by virtue of relying on shared approximation techniques and sources. Note: Benchmark inputs for powr(f) are identical to pow(f). Performance gain over pow on V1 with GCC@15: - SVE powr: 10-12% on subnormal x, 12-13% on x < 0. - SVE powrf: 15% on all x < 0. - AdvSIMD powr: for x < 0, 40% if x subnormal, 60% otherwise. - AdvSIMD powrf: 4% on x subnormals or x < 0.
242 lines
9.2 KiB
C
242 lines
9.2 KiB
C
/* Helper for AdvSIMD single-precision powr
|
|
|
|
Copyright (C) 2025-2026 Free Software Foundation, Inc.
|
|
This file is part of the GNU C Library.
|
|
|
|
The GNU C Library is free software; you can redistribute it and/or
|
|
modify it under the terms of the GNU Lesser General Public
|
|
License as published by the Free Software Foundation; either
|
|
version 2.1 of the License, or (at your option) any later version.
|
|
|
|
The GNU C Library is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Lesser General Public License for more details.
|
|
|
|
You should have received a copy of the GNU Lesser General Public
|
|
License along with the GNU C Library; if not, see
|
|
<https://www.gnu.org/licenses/>. */
|
|
|
|
#include "powf_common.h"
|
|
|
|
#define Log2IdxMask (V_POWF_LOG2_N - 1)
|
|
#define Exp2IdxMask (V_POWF_EXP2_N - 1)
|
|
#define Scale ((double) V_POWF_EXP2_N)
|
|
#define SignBias (1 << (V_POWF_EXP2_TABLE_BITS + 11))
|
|
|
|
#define MantissaMask 0x007fffff
|
|
|
|
static const struct data
|
|
{
|
|
uint32x4_t one, special_bound, sign_bias;
|
|
float32x4_t norm;
|
|
uint32x4_t subnormal_bias;
|
|
uint32x4_t off;
|
|
float32x4_t uflow_bound, oflow_bound;
|
|
uint32x4_t inf;
|
|
float32x4_t nan;
|
|
struct
|
|
{
|
|
double invc, logc;
|
|
} log2_tab[V_POWF_LOG2_N];
|
|
float64x2_t log2_poly[4];
|
|
uint64_t exp2_tab[V_POWF_EXP2_N];
|
|
float64x2_t exp2_poly[3];
|
|
} data = {
|
|
/* Table and polynomial for log2 approximation. */
|
|
.log2_tab = {
|
|
{0x1.6489890582816p+0, -0x1.e960f97b22702p-2 * Scale},
|
|
{0x1.5cf19b35e3472p+0, -0x1.c993406cd4db6p-2 * Scale},
|
|
{0x1.55aac0e956d65p+0, -0x1.aa711d9a7d0f3p-2 * Scale},
|
|
{0x1.4eb0022977e01p+0, -0x1.8bf37bacdce9bp-2 * Scale},
|
|
{0x1.47fcccda1dd1fp+0, -0x1.6e13b3519946ep-2 * Scale},
|
|
{0x1.418ceabab68c1p+0, -0x1.50cb8281e4089p-2 * Scale},
|
|
{0x1.3b5c788f1edb3p+0, -0x1.341504a237e2bp-2 * Scale},
|
|
{0x1.3567de48e9c9ap+0, -0x1.17eaab624ffbbp-2 * Scale},
|
|
{0x1.2fabc80fd19bap+0, -0x1.f88e708f8c853p-3 * Scale},
|
|
{0x1.2a25200ce536bp+0, -0x1.c24b6da113914p-3 * Scale},
|
|
{0x1.24d108e0152e3p+0, -0x1.8d02ee397cb1dp-3 * Scale},
|
|
{0x1.1facd8ab2fbe1p+0, -0x1.58ac1223408b3p-3 * Scale},
|
|
{0x1.1ab614a03efdfp+0, -0x1.253e6fd190e89p-3 * Scale},
|
|
{0x1.15ea6d03af9ffp+0, -0x1.e5641882c12ffp-4 * Scale},
|
|
{0x1.1147b994bb776p+0, -0x1.81fea712926f7p-4 * Scale},
|
|
{0x1.0ccbf650593aap+0, -0x1.203e240de64a3p-4 * Scale},
|
|
{0x1.0875408477302p+0, -0x1.8029b86a78281p-5 * Scale},
|
|
{0x1.0441d42a93328p+0, -0x1.85d713190fb9p-6 * Scale},
|
|
{0x1p+0, 0x0p+0 * Scale},
|
|
{0x1.f1d006c855e86p-1, 0x1.4c1cc07312997p-5 * Scale},
|
|
{0x1.e28c3341aa301p-1, 0x1.5e1848ccec948p-4 * Scale},
|
|
{0x1.d4bdf9aa64747p-1, 0x1.04cfcb7f1196fp-3 * Scale},
|
|
{0x1.c7b45a24e5803p-1, 0x1.582813d463c21p-3 * Scale},
|
|
{0x1.bb5f5eb2ed60ap-1, 0x1.a936fa68760ccp-3 * Scale},
|
|
{0x1.afb0bff8fe6b4p-1, 0x1.f81bc31d6cc4ep-3 * Scale},
|
|
{0x1.a49badf7ab1f5p-1, 0x1.2279a09fae6b1p-2 * Scale},
|
|
{0x1.9a14a111fc4c9p-1, 0x1.47ec0b6df5526p-2 * Scale},
|
|
{0x1.901131f5b2fdcp-1, 0x1.6c71762280f1p-2 * Scale},
|
|
{0x1.8687f73f6d865p-1, 0x1.90155070798dap-2 * Scale},
|
|
{0x1.7d7067eb77986p-1, 0x1.b2e23b1d3068cp-2 * Scale},
|
|
{0x1.74c2c1cf97b65p-1, 0x1.d4e21b0daa86ap-2 * Scale},
|
|
{0x1.6c77f37cff2a1p-1, 0x1.f61e2a2f67f3fp-2 * Scale},
|
|
},
|
|
.log2_poly = { /* rel err: 1.5 * 2^-30. */
|
|
V2 (-0x1.6ff5daa3b3d7cp-2 * Scale),
|
|
V2 (0x1.ec81d03c01aebp-2 * Scale),
|
|
V2 (-0x1.71547bb43f101p-1 * Scale),
|
|
V2 (0x1.7154764a815cbp0 * Scale)
|
|
},
|
|
/* Table and polynomial for exp2 approximation. */
|
|
.exp2_tab = {
|
|
0x3ff0000000000000, 0x3fefd9b0d3158574, 0x3fefb5586cf9890f,
|
|
0x3fef9301d0125b51, 0x3fef72b83c7d517b, 0x3fef54873168b9aa,
|
|
0x3fef387a6e756238, 0x3fef1e9df51fdee1, 0x3fef06fe0a31b715,
|
|
0x3feef1a7373aa9cb, 0x3feedea64c123422, 0x3feece086061892d,
|
|
0x3feebfdad5362a27, 0x3feeb42b569d4f82, 0x3feeab07dd485429,
|
|
0x3feea47eb03a5585, 0x3feea09e667f3bcd, 0x3fee9f75e8ec5f74,
|
|
0x3feea11473eb0187, 0x3feea589994cce13, 0x3feeace5422aa0db,
|
|
0x3feeb737b0cdc5e5, 0x3feec49182a3f090, 0x3feed503b23e255d,
|
|
0x3feee89f995ad3ad, 0x3feeff76f2fb5e47, 0x3fef199bdd85529c,
|
|
0x3fef3720dcef9069, 0x3fef5818dcfba487, 0x3fef7c97337b9b5f,
|
|
0x3fefa4afa2a490da, 0x3fefd0765b6e4540,
|
|
},
|
|
.exp2_poly = { /* rel err: 1.69 * 2^-34. */
|
|
V2 (0x1.c6af84b912394p-5 / Scale / Scale / Scale),
|
|
V2 (0x1.ebfce50fac4f3p-3 / Scale / Scale),
|
|
V2 (0x1.62e42ff0c52d6p-1 / Scale),
|
|
},
|
|
.one = V4 (1),
|
|
.special_bound = V4 (2u * 0x7f800000 - 1),
|
|
.norm = V4 (0x1p23f),
|
|
.subnormal_bias = V4 (0x0b800000), /* 23 << 23. */
|
|
.off = V4 (0x3f35d000),
|
|
.sign_bias = V4 (SignBias),
|
|
.inf = V4 (0x7f800000),
|
|
.nan = V4 (__builtin_nanf ("")),
|
|
/* 2.6 ulp ~ 0.5 + 2^24 (128*Ln2*relerr_log2 + relerr_exp2). */
|
|
.uflow_bound = V4 (-0x1.2cp+12f), /* -150.0 * V_POWF_EXP2_N. */
|
|
.oflow_bound = V4 (0x1p+12f), /* 128.0 * V_POWF_EXP2_N. */
|
|
};
|
|
|
|
/* Check if zero, inf or nan. */
|
|
static inline uint32x4_t
|
|
v_zeroinfnan (const struct data *d, uint32x4_t i)
|
|
{
|
|
return vcgeq_u32 (vsubq_u32 (vaddq_u32 (i, i), d->one), d->special_bound);
|
|
}
|
|
|
|
static inline float64x2_t
|
|
ylogx_core (const struct data *d, float64x2_t iz, float64x2_t k,
|
|
float64x2_t invc, float64x2_t logc, float64x2_t y)
|
|
{
|
|
|
|
/* log2(x) = log1p(z/c-1)/ln2 + log2(c) + k. */
|
|
float64x2_t r = vfmaq_f64 (v_f64 (-1.0), iz, invc);
|
|
float64x2_t y0 = vaddq_f64 (logc, k);
|
|
|
|
/* Polynomial to approximate log1p(r)/ln2. */
|
|
float64x2_t logx = vfmaq_f64 (d->log2_poly[1], r, d->log2_poly[0]);
|
|
logx = vfmaq_f64 (d->log2_poly[2], logx, r);
|
|
logx = vfmaq_f64 (d->log2_poly[3], logx, r);
|
|
logx = vfmaq_f64 (y0, logx, r);
|
|
|
|
return vmulq_f64 (logx, y);
|
|
}
|
|
|
|
static inline float64x2_t
|
|
log2_lookup (const struct data *d, uint32_t i)
|
|
{
|
|
return vld1q_f64 (
|
|
&d->log2_tab[(i >> (23 - V_POWF_LOG2_TABLE_BITS)) & Log2IdxMask].invc);
|
|
}
|
|
|
|
static inline uint64x1_t
|
|
exp2_lookup (const struct data *d, uint64_t i)
|
|
{
|
|
return vld1_u64 (&d->exp2_tab[i & Exp2IdxMask]);
|
|
}
|
|
|
|
static inline float64x2_t
|
|
exp2_core (const struct data *d, float64x2_t ylogx)
|
|
{
|
|
/* N*x = k + r with r in [-1/2, 1/2]. */
|
|
float64x2_t kd = vrndnq_f64 (ylogx);
|
|
int64x2_t ki = vcvtaq_s64_f64 (ylogx);
|
|
float64x2_t r = vsubq_f64 (ylogx, kd);
|
|
|
|
/* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1). */
|
|
uint64x2_t t = vcombine_u64 (exp2_lookup (d, vgetq_lane_s64 (ki, 0)),
|
|
exp2_lookup (d, vgetq_lane_s64 (ki, 1)));
|
|
t = vaddq_u64 (t, vreinterpretq_u64_s64 (
|
|
vshlq_n_s64 (ki, 52 - V_POWF_EXP2_TABLE_BITS)));
|
|
float64x2_t s = vreinterpretq_f64_u64 (t);
|
|
float64x2_t p = vfmaq_f64 (d->exp2_poly[1], r, d->exp2_poly[0]);
|
|
p = vfmaq_f64 (d->exp2_poly[2], r, p);
|
|
p = vfmaq_f64 (s, p, vmulq_f64 (s, r));
|
|
return p;
|
|
}
|
|
|
|
static inline float32x4_t
|
|
powrf_core (const struct data *d, float32x4_t *ylogx, uint32x4_t tmp,
|
|
float32x4_t iz, float32x4_t y, int32x4_t k)
|
|
{
|
|
/* Use double precision for each lane: split input vectors into lo and hi
|
|
halves and promote. */
|
|
float64x2_t tab0 = log2_lookup (d, vgetq_lane_u32 (tmp, 0)),
|
|
tab1 = log2_lookup (d, vgetq_lane_u32 (tmp, 1)),
|
|
tab2 = log2_lookup (d, vgetq_lane_u32 (tmp, 2)),
|
|
tab3 = log2_lookup (d, vgetq_lane_u32 (tmp, 3));
|
|
|
|
float64x2_t iz_lo = vcvt_f64_f32 (vget_low_f32 (iz)),
|
|
iz_hi = vcvt_high_f64_f32 (iz);
|
|
|
|
float64x2_t k_lo = vcvtq_f64_s64 (vmovl_s32 (vget_low_s32 (k))),
|
|
k_hi = vcvtq_f64_s64 (vmovl_high_s32 (k));
|
|
|
|
float64x2_t invc_lo = vzip1q_f64 (tab0, tab1),
|
|
invc_hi = vzip1q_f64 (tab2, tab3),
|
|
logc_lo = vzip2q_f64 (tab0, tab1),
|
|
logc_hi = vzip2q_f64 (tab2, tab3);
|
|
|
|
float64x2_t y_lo = vcvt_f64_f32 (vget_low_f32 (y)),
|
|
y_hi = vcvt_high_f64_f32 (y);
|
|
|
|
float64x2_t ylogx_lo = ylogx_core (d, iz_lo, k_lo, invc_lo, logc_lo, y_lo);
|
|
float64x2_t ylogx_hi = ylogx_core (d, iz_hi, k_hi, invc_hi, logc_hi, y_hi);
|
|
|
|
float32x2_t p_lo = vcvt_f32_f64 (exp2_core (d, ylogx_lo));
|
|
float32x2_t p_hi = vcvt_f32_f64 (exp2_core (d, ylogx_hi));
|
|
|
|
*ylogx = vcombine_f32 (vcvt_f32_f64 (ylogx_lo), vcvt_f32_f64 (ylogx_hi));
|
|
|
|
return vcombine_f32 (p_lo, p_hi);
|
|
}
|
|
|
|
/* Power implementation without assumptions on x or y.
|
|
Evaluate powr(|x|) = exp (y * log(|x|)), and handle
|
|
sign of x using the sign bias.
|
|
Handle underflow and overflow in exponential. */
|
|
static inline float32x4_t
|
|
v_powrf_core (float32x4_t x, float32x4_t y, const struct data *d)
|
|
{
|
|
uint32x4_t ix = vreinterpretq_u32_f32 (x);
|
|
|
|
/* Part of core computation carried in working precision. */
|
|
uint32x4_t tmp = vsubq_u32 (ix, d->off);
|
|
uint32x4_t top = vbicq_u32 (tmp, v_u32 (MantissaMask));
|
|
float32x4_t iz = vreinterpretq_f32_u32 (vsubq_u32 (ix, top));
|
|
int32x4_t k
|
|
= vshrq_n_s32 (vreinterpretq_s32_u32 (top),
|
|
23 - V_POWF_EXP2_TABLE_BITS); /* arithmetic shift. */
|
|
|
|
/* Compute core in extended precision and return intermediate ylogx results
|
|
to handle cases of underflow and overflow in exp. */
|
|
float32x4_t ylogx;
|
|
float32x4_t ret = powrf_core (d, &ylogx, tmp, iz, y, k);
|
|
|
|
/* Handle exp special cases of underflow and overflow. */
|
|
float32x4_t ret_oflow = vreinterpretq_f32_u32 (d->inf);
|
|
float32x4_t ret_uflow = v_f32 (0);
|
|
ret = vbslq_f32 (vcleq_f32 (ylogx, d->uflow_bound), ret_uflow, ret);
|
|
ret = vbslq_f32 (vcgtq_f32 (ylogx, d->oflow_bound), ret_oflow, ret);
|
|
return ret;
|
|
}
|