P3.4: transcribe plugin's own vectorized ln(x) (vlog.cpp, minimax ln(1+x) poly + ln2 range-reduction); correct prior misread — dispatch reaches runtime ln, not FFT; vlog_check ALL OK (rel 2.4e-7)
This commit is contained in:
@@ -19,6 +19,7 @@ add_library(soothe2_dsp SHARED
|
||||
levelpath.cpp
|
||||
phase_table.cpp
|
||||
fftconv.cpp
|
||||
vlog.cpp
|
||||
leveltrack.cpp
|
||||
)
|
||||
|
||||
@@ -26,11 +27,13 @@ add_executable(soothe2_harness harness.cpp)
|
||||
add_executable(twin_check twin_check.cpp)
|
||||
add_executable(tables_check tables_check.cpp)
|
||||
add_executable(fftconv_check fftconv_check.cpp)
|
||||
add_executable(vlog_check vlog_check.cpp)
|
||||
add_executable(leveltrack_check leveltrack_check.cpp)
|
||||
add_executable(levelpath_check levelpath_check.cpp)
|
||||
target_link_libraries(twin_check soothe2_dsp)
|
||||
target_link_libraries(tables_check soothe2_dsp)
|
||||
target_link_libraries(fftconv_check soothe2_dsp)
|
||||
target_link_libraries(vlog_check soothe2_dsp)
|
||||
target_link_libraries(leveltrack_check soothe2_dsp)
|
||||
target_link_libraries(levelpath_check soothe2_dsp)
|
||||
target_link_libraries(soothe2_harness soothe2_dsp)
|
||||
|
||||
@@ -0,0 +1,82 @@
|
||||
#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
|
||||
@@ -0,0 +1,20 @@
|
||||
#pragma once
|
||||
#include <cstdint>
|
||||
|
||||
namespace vlog {
|
||||
|
||||
// Scalar transcription of the vectorized natural-logarithm kernel at 0x1802a24c0
|
||||
// (soothe2 VST3, MSVC x86-64 AVX2, extracted from /tmp/snap_rt.bin).
|
||||
//
|
||||
// The kernel computes dst[i] = ln(src[i]) elementwise in single precision,
|
||||
// using its own minimax polynomial + ln(2) range reduction driven by float
|
||||
// bit-manipulation (the "2/3" magic exponent/mantissa split). A scalar
|
||||
// double-precision Cody-Waite ln handles the slow path for special values.
|
||||
//
|
||||
// NOTE: despite the file name, this kernel is a natural logarithm, NOT an FFT
|
||||
// butterfly. It sits in the spectral (log-magnitude) processing path, not the
|
||||
// complex FFT transform. The signature below mirrors the actual code: float,
|
||||
// out-of-place, real (src != dst is allowed).
|
||||
void log_f32(const float* src, float* dst, uint32_t n);
|
||||
|
||||
} // namespace vlog
|
||||
@@ -0,0 +1,63 @@
|
||||
#include "vlog.hpp"
|
||||
#include <cmath>
|
||||
#include <cstdio>
|
||||
#include <cstdint>
|
||||
#include <vector>
|
||||
|
||||
// Checks the 0x1802a24c0 transcription against the reference natural logarithm.
|
||||
//
|
||||
// The kernel is a vectorized single-precision ln(x), so the reference is
|
||||
// std::log (double) evaluated on the same float inputs. The reported metric is
|
||||
// the max relative error over well-conditioned points (|ln(x)| > 1e-6); for
|
||||
// inputs where ln(x) ~ 0 (x ~ 1) an absolute error is reported instead.
|
||||
// A float minimax polynomial reaches ~1-2 ulp, so a 1e-6 relative gate is the
|
||||
// right tolerance (a double FFT-style 1e-9 gate would be unreachable for float).
|
||||
|
||||
static double rel_err(double a, double b) {
|
||||
double denom = std::abs(b) > 1e-6 ? std::abs(b) : 1.0;
|
||||
return std::abs(a - b) / denom;
|
||||
}
|
||||
|
||||
int main() {
|
||||
const uint32_t N = 8192;
|
||||
std::vector<float> src(N), dst(N);
|
||||
|
||||
// Sweep a wide dynamic range plus a fine neighbourhood around 1.0.
|
||||
uint32_t i = 0;
|
||||
for (int k = -30; k <= 30 && i < N; ++k) {
|
||||
double v = std::pow(10.0, double(k) * 0.5);
|
||||
src[i++] = static_cast<float>(v);
|
||||
}
|
||||
for (int k = -60; k <= 60 && i < N; ++k) {
|
||||
src[i++] = static_cast<float>(1.0 + double(k) * 1e-3);
|
||||
}
|
||||
for (int k = 0; k < 1024 && i < N; ++k) {
|
||||
src[i++] = static_cast<float>(double(k + 1) / 1024.0);
|
||||
}
|
||||
while (i < N) src[i++] = static_cast<float>(i);
|
||||
|
||||
vlog::log_f32(src.data(), dst.data(), N);
|
||||
|
||||
double max_rel = 0.0, max_abs = 0.0;
|
||||
uint32_t rel_i = 0, abs_i = 0;
|
||||
for (uint32_t j = 0; j < N; ++j) {
|
||||
double ref = std::log(static_cast<double>(src[j]));
|
||||
double mine = static_cast<double>(dst[j]);
|
||||
double e = std::abs(mine - ref);
|
||||
if (e > max_abs) { max_abs = e; abs_i = j; }
|
||||
double r = rel_err(mine, ref);
|
||||
if (r > max_rel) { max_rel = r; rel_i = j; }
|
||||
}
|
||||
|
||||
std::printf("n = %u\n", N);
|
||||
std::printf("max relative error = %.6e (at src=%.9g, got %.12g, ref %.12g)\n",
|
||||
max_rel, src[rel_i], dst[rel_i], std::log(static_cast<double>(src[rel_i])));
|
||||
std::printf("max absolute error = %.6e (at src=%.9g)\n", max_abs, src[abs_i]);
|
||||
|
||||
if (max_rel < 1e-6) {
|
||||
std::printf("ALL OK\n");
|
||||
return 0;
|
||||
}
|
||||
std::printf("FAILED\n");
|
||||
return 1;
|
||||
}
|
||||
+22
-15
@@ -95,22 +95,29 @@ These 6 are the final split-radix FFT kernels (stack frames 0x328/0x350/0x7c8/0x
|
||||
vldmxcsr, round-to-nearest 0x1f80) AND x87 control word (fnstcw/fldcw) before the FP loop,
|
||||
then restores. So the transform runs under an explicitly-forced rounding mode.
|
||||
|
||||
### STATUS / scope
|
||||
Full bit-exact FFT = plan generator (0x2f980 split-radix index algebra) + runtime dispatch
|
||||
(3-level tables) + 6 final kernels (each 0x800-0x1200 bytes of vectorized split-radix butterfly)
|
||||
+ twiddle. This is the P3 2-4 week body. Architecture fully mapped; transcription not yet done.
|
||||
### CORRECTION (2026-08-20b): the "final kernels" are vectorized ln(x), NOT FFT
|
||||
The 6 addresses 0x1802a24c0..0x1802ce4a0 are the plugin's OWN vectorized **natural log**,
|
||||
not split-radix FFT butterflies (verified by subagent numeric simulation: matches std::log
|
||||
to float precision). The dispatch chain 0x535a70→…→[0x1826181d8] reaches the plugin's
|
||||
runtime math-function table, not the FFT. Three ln variants live around 0x1802a24c0:
|
||||
- 0x1802a24c0 = AVX2 float ln (minimax poly, range-reduction via 2/3 magic 0x3f2aaaab, ln2)
|
||||
- 0x1802a2fc0 = scalar double ln (Cody-Waite table 128×3 + Taylor), slow path
|
||||
- 0x1802a3260 = second ln variant (9-term, split ln2 hi/lo)
|
||||
The earlier "own vectorized sin/cos" reading was WRONG — the constants (0.333366, −0.250047,
|
||||
…, 0.693147=ln2, 2/3, 0.75) are the minimax coefficients of ln(1+x), not sin/cos.
|
||||
Transcribed to `dsp/vlog.{hpp,cpp}` (namespace vlog, `log_f32(src,dst,n)`), check ALL OK
|
||||
(max rel err 2.4e-7 vs std::log). This is the exact log used by level-path (logf·8.6859
|
||||
for dB, log/gamma in LUT curve), so it feeds bit-exact level→dB.
|
||||
|
||||
### Final kernel 0x1802a24c0 — own vectorized sin/cos, NOT twiddle table
|
||||
The big-N kernels compute trig ON THE FLY via an AVX polynomial (vfmadd231ps/213ps chain) with
|
||||
range reduction done by float-bit tricks (vpaddd/vpsubd/vpsrad $0x17 = 23-bit shift = exponent
|
||||
extract, vpcmpgtd, vblendvps). Polynomial constants (rodata, float, broadcast x8):
|
||||
0x181f81fa0 = 0.333366 0x181f81fc0 = -0.250047 0x181f81fe0 = 0.198225
|
||||
0x181f82000 = -0.164625 0x181f82020 = 0.169649 0x181f82040 = -0.151772
|
||||
0x181f82100 = 0.666667 0x181f821c0 = 0.693147 (= ln 2) 0x181f822c0 = 0.75
|
||||
(0.693147 = ln 2 ⇒ exp-based sin/cos, SVML-like vector math library.)
|
||||
=> Small-N path uses the twiddle sin-table (0x39b00); big-N path uses this own polynomial.
|
||||
Transcribing this bit-exactly = reimplementing a vector math sin/cos (Cody-Waite + poly + exp),
|
||||
plus the split-radix butterfly + integer reorder — the multi-week P3 body.
|
||||
### REAL FFT pieces (still to transcribe):
|
||||
- stage kernels FUN_18000bfc0 / 18000c5e0 (butterfly) — see above
|
||||
- plan generator FUN_18002f980 (split-radix index algebra)
|
||||
- actual twiddle usage for small N via FUN_180039b00 sin-table
|
||||
So P3 scope is unchanged (butterfly + plan gen), NOT the ln kernels.
|
||||
|
||||
## STATUS / scope
|
||||
Full bit-exact FFT = plan generator (0x2f980 split-radix index algebra) + stage kernels
|
||||
(0xbfc0/0xc5e0 butterfly) + twiddle. ln (vlog) now done. FFT butterfly + plan gen remain.
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user