diff --git a/dsp/CMakeLists.txt b/dsp/CMakeLists.txt index 334687a..a80bdb2 100644 --- a/dsp/CMakeLists.txt +++ b/dsp/CMakeLists.txt @@ -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) diff --git a/dsp/vlog.cpp b/dsp/vlog.cpp new file mode 100644 index 0000000..bd9210e --- /dev/null +++ b/dsp/vlog.cpp @@ -0,0 +1,82 @@ +#include "vlog.hpp" +#include +#include + +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(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(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(std::log(static_cast(y))); + } + } +} + +} // namespace vlog diff --git a/dsp/vlog.hpp b/dsp/vlog.hpp new file mode 100644 index 0000000..18a1f98 --- /dev/null +++ b/dsp/vlog.hpp @@ -0,0 +1,20 @@ +#pragma once +#include + +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 diff --git a/dsp/vlog_check.cpp b/dsp/vlog_check.cpp new file mode 100644 index 0000000..7776390 --- /dev/null +++ b/dsp/vlog_check.cpp @@ -0,0 +1,63 @@ +#include "vlog.hpp" +#include +#include +#include +#include + +// 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 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(v); + } + for (int k = -60; k <= 60 && i < N; ++k) { + src[i++] = static_cast(1.0 + double(k) * 1e-3); + } + for (int k = 0; k < 1024 && i < N; ++k) { + src[i++] = static_cast(double(k + 1) / 1024.0); + } + while (i < N) src[i++] = static_cast(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(src[j])); + double mine = static_cast(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(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; +} diff --git a/notes_giant_fft.md b/notes_giant_fft.md index bcd1bdb..408b23a 100644 --- a/notes_giant_fft.md +++ b/notes_giant_fft.md @@ -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.