83 lines
3.3 KiB
C++
83 lines
3.3 KiB
C++
#include "vlog.hpp"
|
|
#include <cmath>
|
|
#include <cstring>
|
|
|
|
namespace vlog {
|
|
|
|
namespace {
|
|
|
|
// Float bit patterns used by the AVX kernel's integer lane ops.
|
|
constexpr uint32_t C_2_OVER_3_BITS = 0x3f2aaaab; // bits of 2/3f (0.6666666865f)
|
|
constexpr uint32_t C_MANTISSA_MASK = 0x007fffff; // 2^23 - 1 (low 23 mantissa bits)
|
|
|
|
// ln(2) as a single float constant (0x3f317218).
|
|
constexpr float C_LN2 = 0.6931471824645996f;
|
|
|
|
// Minimax polynomial coefficients for ln(1+x) on the reduced interval
|
|
// x in [-1/3, 1/3), read straight from the binary:
|
|
// ln(1+x) ~= x + x^2 * (c7 + c6*x + c5*x^2 + c4*x^3 + c3*x^4 + c2*x^5 + c1*x^6)
|
|
constexpr float C_P1 = -0.15177205204963684f; // 0x181f82040
|
|
constexpr float C_P2 = 0.16964881122112274f; // 0x181f82020
|
|
constexpr float C_P3 = -0.16462457180023193f; // 0x181f82000
|
|
constexpr float C_P4 = 0.19822503626346588f; // 0x181f81fe0
|
|
constexpr float C_P5 = -0.25004664063453674f; // 0x181f81fc0
|
|
constexpr float C_P6 = 0.33336564898490906f; // 0x181f81fa0
|
|
constexpr float C_P7 = -0.5f; // 0x181f81f80
|
|
|
|
inline float fma_f(float a, float b, float c) { return __builtin_fmaf(a, b, c); }
|
|
|
|
inline uint32_t bit_u32(float f) { uint32_t u; std::memcpy(&u, &f, sizeof u); return u; }
|
|
inline float bit_f32(uint32_t u) { float f; std::memcpy(&f, &u, sizeof f); return f; }
|
|
|
|
// Fast path: input is a positive normal float -> ln(x).
|
|
//
|
|
// Range reduction via the "2/3" magic (equivalent to the kernel's
|
|
// vpsubd/vpsrad 0x17/vpand/vpaddd): the mantissa is folded into
|
|
// m in [2/3, 4/3) and the integer exponent e is recovered, so that
|
|
// y = m * 2^e => ln(y) = e*ln2 + ln(m),
|
|
// with ln(m) = ln(1+x), x = m-1 in [-1/3, 1/3), from the minimax polynomial.
|
|
inline float ln_fast(float y) {
|
|
uint32_t b = bit_u32(y);
|
|
uint32_t t = b - C_2_OVER_3_BITS; // vpsubd (wrapping)
|
|
int32_t e = static_cast<int32_t>(t) >> 23; // vpsrad 0x17 (exponent)
|
|
uint32_t mb = (t & C_MANTISSA_MASK) + C_2_OVER_3_BITS; // vpand + vpaddd
|
|
float m = bit_f32(mb);
|
|
float x = m - 1.0f; // vsubps (x in [-1/3, 1/3))
|
|
float e_ = static_cast<float>(e); // vcvtdq2ps
|
|
|
|
// Horner evaluation (mirrors the vfmadd231ps/vfmadd213ps chain).
|
|
float p = C_P2;
|
|
p = fma_f(x, C_P1, p);
|
|
p = fma_f(p, x, C_P3);
|
|
p = fma_f(p, x, C_P4);
|
|
p = fma_f(p, x, C_P5);
|
|
p = fma_f(p, x, C_P6);
|
|
p = fma_f(p, x, C_P7);
|
|
|
|
float q = x * p; // x*P(x)
|
|
q = fma_f(q, x, x); // x^2*P(x) + x ~= ln(1+x) = ln(m)
|
|
return fma_f(e_, C_LN2, q); // e*ln2 + ln(m) == ln(y)
|
|
}
|
|
|
|
} // namespace
|
|
|
|
void log_f32(const float* src, float* dst, uint32_t n) {
|
|
for (uint32_t i = 0; i < n; ++i) {
|
|
float y = src[i];
|
|
uint32_t b = bit_u32(y);
|
|
|
|
// Fast path iff the lane is a positive normal float:
|
|
// (int32)(b + 0x00800000) >= 0x01000000 <=> b in [0x00800000, 0x7f7fffff].
|
|
if (b >= 0x00800000u && b <= 0x7f7fffffu) {
|
|
dst[i] = ln_fast(y);
|
|
} else {
|
|
// Slow path: zero/denormal/negative/Inf/NaN. The kernel dispatches to
|
|
// a scalar double-precision Cody-Waite ln (0x1802a2fc0); std::log is the
|
|
// structurally equivalent reference for these edge inputs.
|
|
dst[i] = static_cast<float>(std::log(static_cast<double>(y)));
|
|
}
|
|
}
|
|
}
|
|
|
|
} // namespace vlog
|