From 51c46db0fe9793083faa9b03f355cdc0e80e8df4 Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Thu, 17 Sep 2026 18:20:11 +0000 Subject: [PATCH 1/4] Add the Stage 2a FFT test suite: Ooura parity gate, independent oracle, RT guard Written first, against the vendored C, so the wave-2 port works test-first (docs/audit-fft-and-code-smells.md, Part 3 Stage 2a and Part 9). - tests/test_fft_parity_ooura.cpp + target tap_dsp_fft_parity: memcmp bit identity of basic_real_fft against raw rdft/rdft_f from a private reference build of fftsg.c/fftsg_float.c, both sides compiled with -ffp-contract=off (MSVC: nothing), forward and inverse, N = 4..65536 plus 2^20, on five materials. ooura_ref and engine_under_test are each re-pointable in one place. TAP_DSP_PARITY_MAX_N (default 2^20) caps the sizes for emulated targets. A second target, tap_dsp_fft_parity_default_flags, prints the measured max-ulp deviation per N at default flags and never fails (Part 6, N3). - tests/test_fft_oracle.cpp: closed-form vectors (impulse, DC, Nyquist, on-bin cosine/sine with the +i convention, two-tone) and their unnormalized inverses, to a tolerance derived from Higham's FFT bound (4 * eps * log2 N * ||y||_2, measured 0.25/0.17 of it for double/float); plus a double-double (TwoSum/TwoProd) compensated DFT for N <= 256, with a profile extension point for the fixed-point stage. - tests/test_fft_rt.cpp: static_assert(noexcept) on the four transforms and size queries; a global operator new/delete counting guard proving the four calls allocate nothing at N = 512 and 4096, on the first call after construction and in steady state; copy and copy-assignment bit-identical. - tests/support/signals.h: one xorshift32, random_signal, tone (exact modular phase reduction) and db_from_ratio for the new files. - tests/CMakeLists.txt: additive, delimited block at the end. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- tests/CMakeLists.txt | 69 ++++ tests/support/signals.h | 83 ++++ tests/test_fft_oracle.cpp | 667 ++++++++++++++++++++++++++++++++ tests/test_fft_parity_ooura.cpp | 398 +++++++++++++++++++ tests/test_fft_rt.cpp | 386 ++++++++++++++++++ 5 files changed, 1603 insertions(+) create mode 100644 tests/support/signals.h create mode 100644 tests/test_fft_oracle.cpp create mode 100644 tests/test_fft_parity_ooura.cpp create mode 100644 tests/test_fft_rt.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index eaaa14b..4872cb7 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -140,3 +140,72 @@ tap_dsp_add_gtest_executable(tap_dsp_tests target_link_libraries(tap_dsp_tests PRIVATE tap::dsp tap_dsp_warnings) + +# ============================================================================== +# BEGIN Stage 2a additions (docs/audit-fft-and-code-smells.md, Part 9): the +# independent oracle and the real-time guard join the main executable; the +# Ooura bit-identity gate gets its own target so it can own its compiler flags. +# Keep this block self-contained and at the end of the file. +# ============================================================================== + +target_sources(tap_dsp_tests PRIVATE + test_fft_oracle.cpp + test_fft_rt.cpp) + +# ------------------------------------------------------------------------------ +# The reference copy of the vendored C, compiled with fp-contraction OFF. Bit +# identity between the C and its C++ transliteration holds only if both are +# compiled to the same sequence of IEEE operations, and CMake gives neither +# side an fp-contract setting: gcc contracts C in gnu17 mode, g++ contracts +# C++ in every mode, clang contracts everywhere (Part 4, "fp-contraction +# policy"). MSVC's default /fp:precise does not contract, so nothing is passed +# there. This library is private to the tests and is not the tap_dsp_fft the +# library links; the parity executable links THIS and not tap::dsp, so its +# rdft/rdft_f resolve here and no backend define (vDSP on Apple) reaches it. +# ------------------------------------------------------------------------------ +set(_tap_dsp_nocontract_c "$<$:-ffp-contract=off>") +set(_tap_dsp_nocontract_cxx "$<$:-ffp-contract=off>") + +add_library(tap_dsp_fft_reference_nocontract STATIC + ${PROJECT_SOURCE_DIR}/third_party/ooura/fftsg.c + ${PROJECT_SOURCE_DIR}/third_party/ooura/fftsg_float.c) +target_compile_options(tap_dsp_fft_reference_nocontract PRIVATE ${_tap_dsp_nocontract_c}) + +# Largest transform the parity suite runs. The default covers every consumer +# geometry (2^20); emulated targets set 4096 (Part 10) to fit in RAM. +set(TAP_DSP_PARITY_MAX_N 1048576 CACHE STRING + "Largest FFT size the Ooura parity suite runs (power of two; 4096 on emulated targets)") + +# THE GATE: -ffp-contract=off on both sides, memcmp bit identity. +add_executable(tap_dsp_fft_parity test_fft_parity_ooura.cpp) +target_include_directories(tap_dsp_fft_parity PRIVATE ${PROJECT_SOURCE_DIR}/include) +target_compile_features(tap_dsp_fft_parity PRIVATE cxx_std_20) +target_compile_options(tap_dsp_fft_parity PRIVATE ${_tap_dsp_nocontract_cxx}) +target_compile_definitions(tap_dsp_fft_parity PRIVATE TAP_DSP_PARITY_MAX_N=${TAP_DSP_PARITY_MAX_N}) +target_link_libraries(tap_dsp_fft_parity PRIVATE + tap_dsp_fft_reference_nocontract + tap_dsp_warnings + GTest::gtest_main) +gtest_discover_tests(tap_dsp_fft_parity DISCOVERY_TIMEOUT 120 PROPERTIES TIMEOUT 900 LABELS parity) + +# INFORMATIONAL: the same source at default flags on both sides (the C here is +# the library's own tap_dsp_fft build), printing the measured max-ulp deviation +# per N and never failing. It is not a pass/fail gate with an assumed bound +# because a different fusion choice per butterfly stage accumulates over +# log2 N stages; the number goes on the record per platform instead (Part 6, +# item N3). +add_executable(tap_dsp_fft_parity_default_flags test_fft_parity_ooura.cpp) +target_include_directories(tap_dsp_fft_parity_default_flags PRIVATE ${PROJECT_SOURCE_DIR}/include) +target_compile_features(tap_dsp_fft_parity_default_flags PRIVATE cxx_std_20) +target_compile_definitions(tap_dsp_fft_parity_default_flags PRIVATE + TAP_DSP_PARITY_INFORMATIONAL + TAP_DSP_PARITY_MAX_N=${TAP_DSP_PARITY_MAX_N}) +target_link_libraries(tap_dsp_fft_parity_default_flags PRIVATE + tap::dsp_fft + tap_dsp_warnings + GTest::gtest_main) +gtest_discover_tests(tap_dsp_fft_parity_default_flags DISCOVERY_TIMEOUT 120 PROPERTIES TIMEOUT 900 LABELS parity) + +# ============================================================================== +# END Stage 2a additions +# ============================================================================== diff --git a/tests/support/signals.h b/tests/support/signals.h new file mode 100644 index 0000000..6f6963d --- /dev/null +++ b/tests/support/signals.h @@ -0,0 +1,83 @@ +/// @file signals.h +/// @brief Shared test-signal generators for the DspTap test battery. +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// One xorshift32, one random_signal, one tone synthesizer and one dB helper, +// so a new test file does not grow its own copy (three already exist in +// test_fft.cpp, test_fft_backend.cpp and test_nn.cpp; migrating them is the +// Stage 6 hygiene item in docs/audit-fft-and-code-smells.md, not this file's +// job). Everything here is deterministic — fixed seeds, no wall clock, no +// filesystem — so it can run unchanged on the bare-metal QEMU legs. + +#pragma once + +#include +#include +#include +#include +#include + +namespace tap::dsp::test { + + /// Marsaglia's xorshift32 (Journal of Statistical Software 8(14), 2003), + /// the generator every existing copy in this repo already uses. The state + /// must be non-zero; a zero seed is replaced by a fixed non-zero constant. + class xorshift32 { + public: + explicit xorshift32(std::uint32_t seed) noexcept + : m_s(seed != 0u ? seed : 0x9E3779B9u) {} + + /// Next raw 32-bit state. + std::uint32_t next_u32() noexcept { + m_s ^= m_s << 13; + m_s ^= m_s >> 17; + m_s ^= m_s << 5; + return m_s; + } + + /// Uniform in [-1, 1), computed in double so float and double signals + /// drawn from the same seed are the same values up to the final cast. + double next_unit() noexcept { return static_cast(next_u32()) / 2147483648.0 - 1.0; } + + private: + std::uint32_t m_s; + }; + + /// n samples uniform in [-amplitude, amplitude), from a fixed seed. + template + std::vector random_signal(std::size_t n, std::uint32_t seed, double amplitude = 1.0) { + xorshift32 rng(seed); + std::vector x(n); + for (auto& v : x) { + v = static_cast(amplitude * rng.next_unit()); + } + return x; + } + + /// amplitude * cos(2*pi*bin*j/n + phase) for j in [0, n). + /// + /// The angle is formed as 2*pi * fmod(bin*j, n) / n, not as (2*pi*bin/n) * j: + /// with a rounded per-sample increment the phase error grows like j*ulp, + /// and at n = 65536 a supposedly on-bin cosine is ~1e-11 off, a thousand + /// times the engine's own error — which is what a closed-form oracle would + /// then wrongly report. fmod is exact, bin*j is exact for integer bins below + /// 2^53, and the one division by a power-of-two n is exact, so an integer-bin + /// tone here is the sampled cosine to one rounding of cos() per sample. + template + std::vector tone(std::size_t n, double bin, double amplitude, double phase = 0.0) { + std::vector x(n); + const double period = static_cast(n); + for (std::size_t j = 0; j < n; ++j) { + const double turns = std::fmod(bin * static_cast(j), period) / period; + x[j] = static_cast(amplitude * std::cos(2.0 * std::numbers::pi * turns + phase)); + } + return x; + } + + /// 20*log10(ratio) of an amplitude ratio; -inf for zero, as std::log10 gives. + inline double db_from_ratio(double ratio) noexcept { + return 20.0 * std::log10(ratio); + } + +} // namespace tap::dsp::test diff --git a/tests/test_fft_oracle.cpp b/tests/test_fft_oracle.cpp new file mode 100644 index 0000000..754da1e --- /dev/null +++ b/tests/test_fft_oracle.cpp @@ -0,0 +1,667 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// THE INDEPENDENT ORACLE for tap::dsp::basic_real_fft (Stage 2a of +// docs/audit-fft-and-code-smells.md; Part 9 names this file). +// +// test_fft_parity_ooura.cpp proves the port IS Ooura, bit for bit. It cannot +// prove Ooura is a Fourier transform: if the reference C and the port shared a +// defect, parity would be green. This file answers that with two references +// that share nothing with fftsg.c: +// +// 1. Closed-form vectors whose transforms are known exactly — an impulse +// (flat spectrum), DC (N in slot 0), the Nyquist alternation (N in slot +// 1), an on-bin cosine and sine at bin k (+N/2 in the real, resp. IMAG, +// part of bin k — the plus sign in the imaginary part is the documented +// W = exp(+2*pi*i/N) convention, the conjugate of the engineering DFT), +// and a two-tone superposition — and their unnormalized inverses, whose +// gain is N/2 per the contract in fft.h. Checked to a tolerance DERIVED +// from the sample type's epsilon and log2 N (see tolerance() below). +// +// 2. A compensated-summation DFT for N <= 256: the O(N^2) definition with +// double-double accumulation (TwoSum / TwoProd, Dekker 1971 and Knuth +// TAOCP vol. 2 4.2.2; Shewchuk 1997 for the error-free transformations) +// and double-double twiddles from a Taylor series after exact octant +// reduction on the integer bin*sample index. Not `long double`: that is +// `double` on MSVC and on Apple arm64 (Part 6, item N19), so it would be +// no oracle at all on two of the three hosted CI legs. +// +// The suite is typed over float and double now. The fixed-point stage +// (Stage 3, Part 7) extends `profile` below with int16_t and int32_t, +// supplying the profile's documented forward scale and the conversion into +// the double domain; the tests themselves are written against that trait so +// nothing else changes. + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "support/signals.h" +#include "tap/dsp/fft.h" + +namespace { + + // ------------------------------------------------------------------------ + // Per-profile traits — THE EXTENSION POINT FOR THE FIXED-POINT STAGE. + // + // A profile says how to get a sample into the double domain, what its + // unit roundoff is, and what scale the engine's forward / unnormalized + // inverse carry relative to the mathematical DFT (1 and 1 for the float + // profiles; the Q15/Q31 profiles document X/N under `fixed` scaling and + // will state their own inverse scale). Stage 3 adds: + // template <> struct profile { ... }; + // template <> struct profile { ... }; + // and appends the two types to oracle_types. The tolerance for those + // profiles is a quantization-noise number from Part 7, not epsilon-based; + // tolerance() dispatches on the trait so that is a local change too. + // ------------------------------------------------------------------------ + template + struct profile; + + template <> + struct profile { + static constexpr double k_epsilon = std::numeric_limits::epsilon(); + static constexpr double k_forward_scale = 1.0; ///< engine forward = k * DFT + static constexpr double k_inverse_scale = 1.0; ///< engine inverse = k * Ooura's unnormalized inverse + static double to_double(double v) { return v; } + static double from_double(double v) { return v; } + }; + + template <> + struct profile { + static constexpr double k_epsilon = std::numeric_limits::epsilon(); + static constexpr double k_forward_scale = 1.0; + static constexpr double k_inverse_scale = 1.0; + static double to_double(float v) { return static_cast(v); } + static float from_double(double v) { return static_cast(v); } + }; + + using oracle_types = ::testing::Types; + + // ------------------------------------------------------------------------ + // Tolerance. + // + // Higham, Accuracy and Stability of Numerical Algorithms (2nd ed., 2002), + // Theorem 24.2: a radix-2 FFT of length N = 2^t whose twiddles carry + // relative error at most mu computes y_hat with + // + // || y_hat - y ||_2 <= t * eta / (1 - t * eta) * || y ||_2, + // eta = mu + gamma_4 (1 + mu), gamma_4 ~ 4u, + // + // where u is the unit roundoff (epsilon / 2). Ooura's tables are computed + // from libm and stored in Sample, so mu <= u, hence eta <= 5u + O(u^2) + // and, for t*eta << 1, the 2-norm error is below 5 * u * log2(N) * ||y||_2. + // The largest single-element error cannot exceed the 2-norm, and a + // measured (not remembered) value for || y ||_2 comes free with the exact + // answer. Two further terms ride on top: casting the closed-form input to + // Sample perturbs each x_j by at most u|x_j|, which by Parseval moves y by + // at most u * ||y||_2 in 2-norm; and Ooura is split-radix, not the radix-2 + // of the theorem, so its constant differs in a factor of order one. The + // bound is therefore taken as k_higham_constant * epsilon * log2(N) * + // ||y||_2 with k_higham_constant = 4 (epsilon = 2u, so 4 * epsilon = 8u: + // the 5u of the theorem, the u of the input cast, and 2u of margin for + // the split-radix constant). + // + // Measured against that bound on x86-64 Linux (GCC 13, -O3, glibc libm, + // Ooura as the engine): the largest ratio |error| / (epsilon * log2(N) * + // ||y||_2) over every test in this file was 0.25 for double and 0.17 for + // float, i.e. the engine sits 16x (double) and 23x (float) inside the + // derived bound of 4. That margin is deliberate: this is an oracle for + // correctness, not a precision ratchet; the measured-number pins live in + // test_fft.cpp and test_fft_backend.cpp. + // ------------------------------------------------------------------------ + constexpr double k_higham_constant = 4.0; + + template + double tolerance(std::size_t n, double spectrum_norm2) { + const double log2n = std::log2(static_cast(n)); + return k_higham_constant * profile::k_epsilon * log2n * spectrum_norm2; + } + + /// ||y||_2 of the full complex spectrum of x, from Parseval: sqrt(N) * ||x||_2. + double spectrum_norm2(const std::vector& x) { + double e = 0.0; + for (const double v : x) { + e += v * v; + } + return std::sqrt(static_cast(x.size()) * e); + } + + template + std::vector to_doubles(const std::vector& x) { + std::vector d(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) { + d[i] = profile::to_double(x[i]); + } + return d; + } + + template + std::vector from_doubles(const std::vector& x) { + std::vector s(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) { + s[i] = profile::from_double(x[i]); + } + return s; + } + + // ------------------------------------------------------------------------ + // Double-double arithmetic: an unevaluated sum hi + lo with |lo| <= ulp(hi)/2, + // ~106 bits of significand. Only what the DFT below needs. + // ------------------------------------------------------------------------ + struct dd { + double hi; + double lo; + }; + + /// Knuth TwoSum: s + e == a + b exactly, no ordering assumption. + inline dd two_sum(double a, double b) { + const double s = a + b; + const double bb = s - a; + const double e = (a - (s - bb)) + (b - bb); + return {s, e}; + } + + /// Dekker FastTwoSum, valid when |a| >= |b|. + inline dd quick_two_sum(double a, double b) { + const double s = a + b; + const double e = b - (s - a); + return {s, e}; + } + + /// TwoProd via a correctly rounded fma: p + e == a * b exactly. std::fma + /// is required to round once, whatever the hardware. + inline dd two_prod(double a, double b) { + const double p = a * b; + const double e = std::fma(a, b, -p); + return {p, e}; + } + + inline dd dd_add(dd a, dd b) { + dd s = two_sum(a.hi, b.hi); + dd t = two_sum(a.lo, b.lo); + s.lo += t.hi; + s = quick_two_sum(s.hi, s.lo); + s.lo += t.lo; + return quick_two_sum(s.hi, s.lo); + } + + inline dd dd_neg(dd a) { + return {-a.hi, -a.lo}; + } + + inline dd dd_mul(dd a, dd b) { + dd p = two_prod(a.hi, b.hi); + p.lo += a.hi * b.lo + a.lo * b.hi; + return quick_two_sum(p.hi, p.lo); + } + + inline dd dd_mul_d(dd a, double b) { + dd p = two_prod(a.hi, b); + p.lo += a.lo * b; + return quick_two_sum(p.hi, p.lo); + } + + inline dd dd_div_d(dd a, double b) { + const double q1 = a.hi / b; + const dd r = dd_add(a, dd_neg(dd_mul_d({b, 0.0}, q1))); + const double q2 = r.hi / b; + return quick_two_sum(q1, q2); + } + + inline double dd_round(dd a) { + return a.hi + a.lo; + } + + // pi as a double-double: the standard 3.141592653589793 + 1.2246467991473532e-16. + constexpr dd k_pi_dd{3.141592653589793, 1.2246467991473532e-16}; + + /// sin and cos of |theta| <= pi/4 by Taylor series in double-double. Terms + /// fall below 1e-34 relative by k = 15 for |theta| <= pi/4; 20 terms are + /// summed for margin, largest first, which double-double handles exactly + /// enough for a 1e-30 result. The oracle needs ~1e-20, so this is ample. + struct dd_sincos { + dd sin; + dd cos; + }; + + dd_sincos sincos_small(dd theta) { + const dd theta2 = dd_mul(theta, theta); + dd s_term = theta; + dd c_term = {1.0, 0.0}; + dd s = s_term; + dd c = c_term; + for (int k = 1; k <= 20; ++k) { + // sin: term_k = -term_{k-1} * theta^2 / ((2k)(2k+1)) + // cos: term_k = -term_{k-1} * theta^2 / ((2k-1)(2k)) + s_term = dd_neg(dd_div_d(dd_mul(s_term, theta2), static_cast((2 * k) * (2 * k + 1)))); + c_term = dd_neg(dd_div_d(dd_mul(c_term, theta2), static_cast((2 * k - 1) * (2 * k)))); + s = dd_add(s, s_term); + c = dd_add(c, c_term); + } + return {s, c}; + } + + /// (cos, sin)(2*pi*m/n) for integer m in [0, n) and n a power of two >= 4, + /// by exact octant reduction on the integers (m and n are exact, m/n is + /// exact because n is a power of two) followed by the small-angle series. + dd_sincos twiddle(std::size_t m, std::size_t n) { + const std::size_t half = n / 2; + const std::size_t quarter = n / 4; + bool negate = false; + bool rotate = false; // (cos, sin) -> (-sin, cos) + m %= n; + if (m >= half) { + m -= half; + negate = true; + } + if (m >= quarter) { + m -= quarter; + rotate = true; + } + // m in [0, quarter): theta in [0, pi/2). Fold to [0, pi/4] via + // cos(theta) = sin(pi/2 - theta). + bool swap = false; + if (2 * m > quarter) { + m = quarter - m; + swap = true; + } + const double fraction = static_cast(m) / static_cast(n); // exact + const dd theta = dd_mul_d({2.0 * k_pi_dd.hi, 2.0 * k_pi_dd.lo}, fraction); + const dd_sincos sc = sincos_small(theta); + dd c = swap ? sc.sin : sc.cos; + dd s = swap ? sc.cos : sc.sin; + if (rotate) { + const dd t = c; + c = dd_neg(s); + s = t; + } + if (negate) { + c = dd_neg(c); + s = dd_neg(s); + } + return {s, c}; + } + + // ------------------------------------------------------------------------ + // The compensated DFT, in Ooura's packing and sign convention (fft.h): + // forward: a[0] = Re X[0], a[1] = Re X[N/2], a[2k] = Re X[k], a[2k+1] = Im X[k], + // X[k] = sum_j x[j] * exp(+2*pi*i*j*k/N) + // inverse (unnormalized, gain N/2 on a round trip): + // x[k] = (R[0] + R[N/2] * (-1)^k) / 2 + // + sum_{j=1}^{N/2-1} (R[j] cos(2*pi*j*k/N) + I[j] sin(2*pi*j*k/N)) + // The inverse formula is fftsg.c's own statement of what rdft(n, -1, ...) + // computes; both are restated from the definition, not from the code path. + // ------------------------------------------------------------------------ + class compensated_dft { + public: + explicit compensated_dft(std::size_t n) + : m_n(n) + , m_twiddles(n) { + for (std::size_t m = 0; m < n; ++m) { + m_twiddles[m] = twiddle(m, n); + } + } + + std::vector forward(const std::vector& x) const { + std::vector a(m_n, 0.0); + for (std::size_t k = 0; k <= m_n / 2; ++k) { + dd re{0.0, 0.0}; + dd im{0.0, 0.0}; + for (std::size_t j = 0; j < m_n; ++j) { + const dd_sincos& w = m_twiddles[(j * k) % m_n]; + re = dd_add(re, dd_mul_d(w.cos, x[j])); + im = dd_add(im, dd_mul_d(w.sin, x[j])); + } + if (k == 0) { + a[0] = dd_round(re); + } + else if (k == m_n / 2) { + a[1] = dd_round(re); + } + else { + a[2 * k] = dd_round(re); + a[2 * k + 1] = dd_round(im); + } + } + return a; + } + + std::vector inverse_unnormalized(const std::vector& a) const { + std::vector x(m_n, 0.0); + for (std::size_t k = 0; k < m_n; ++k) { + const double edge = (k % 2 == 0) ? a[0] + a[1] : a[0] - a[1]; + dd acc = dd_mul_d({edge, 0.0}, 0.5); + for (std::size_t j = 1; j < m_n / 2; ++j) { + const dd_sincos& w = m_twiddles[(j * k) % m_n]; + acc = dd_add(acc, dd_mul_d(w.cos, a[2 * j])); + acc = dd_add(acc, dd_mul_d(w.sin, a[2 * j + 1])); + } + x[k] = dd_round(acc); + } + return x; + } + + const dd_sincos& twiddle_at(std::size_t m) const { return m_twiddles[m % m_n]; } + + private: + std::size_t m_n; + std::vector m_twiddles; + }; + + // ------------------------------------------------------------------------ + // Shared checkers. + // ------------------------------------------------------------------------ + template + void expect_close(const std::vector& got, const std::vector& expected, double tol, const char* what, + std::size_t n) { + ASSERT_EQ(got.size(), expected.size()); + for (std::size_t i = 0; i < got.size(); ++i) { + const double g = profile::to_double(got[i]); + ASSERT_NEAR(g, expected[i], tol) << what << " N=" << n << " index " << i; + } + } + + std::vector closed_form_sizes() { + std::vector s; + for (std::size_t n = 4; n <= 65536; n *= 2) { + s.push_back(n); + } + return s; + } + + std::vector dft_sizes() { + std::vector s; + for (std::size_t n = 4; n <= 256; n *= 2) { + s.push_back(n); + } + return s; + } + + /// Bin used for the on-bin materials: n/8, or 1 below n = 8. Always in + /// [1, n/2), i.e. a genuinely complex bin, never DC or Nyquist. + std::size_t tone_bin(std::size_t n) { + return std::max(1, n / 8); + } + + // Runs the engine forward on x (given in the double domain, cast to the + // profile) and checks it against the exact spectrum `expected`. + template + void check_forward(const std::vector& x, const std::vector& expected, const char* what) { + const std::size_t n = x.size(); + tap::dsp::basic_real_fft fft(n); + std::vector buf = from_doubles(x); + fft.forward_inplace(buf.data()); + std::vector scaled = expected; + for (double& v : scaled) { + v *= profile::k_forward_scale; + } + expect_close(buf, scaled, tolerance(n, spectrum_norm2(x)), what, n); + } + + // Runs the engine's UNNORMALIZED inverse on the packed spectrum `a` and + // checks it against `expected` (which already carries the N/2 gain). + template + void check_inverse(const std::vector& a, const std::vector& expected, const char* what) { + const std::size_t n = a.size(); + tap::dsp::basic_real_fft fft(n); + std::vector buf = from_doubles(a); + fft.inverse_inplace(buf.data()); + std::vector scaled = expected; + for (double& v : scaled) { + v *= profile::k_inverse_scale; + } + // ||expected||_2 plays the role of ||y||_2 for the inverse direction; + // the packed spectrum's 2-norm is within sqrt(2) of the true one, and + // the N/2 gain means ||expected||_2 >= ||a||_2 * sqrt(N)/2 for these + // vectors, so the bound stays of the same form. + double e = 0.0; + for (const double v : expected) { + e += v * v; + } + expect_close(buf, scaled, tolerance(n, std::sqrt(e)), what, n); + } + + template + class fft_oracle_test : public ::testing::Test {}; + TYPED_TEST_SUITE(fft_oracle_test, oracle_types); + + // ======================================================================== + // Closed forms, forward. + // ======================================================================== + + TYPED_TEST(fft_oracle_test, ImpulseHasFlatSpectrum) { + for (const std::size_t n : closed_form_sizes()) { + std::vector x(n, 0.0); + x[0] = 1.0; + std::vector expected(n, 0.0); + expected[0] = 1.0; // DC + expected[1] = 1.0; // Nyquist + for (std::size_t k = 1; k < n / 2; ++k) { + expected[2 * k] = 1.0; + } + check_forward(x, expected, "impulse"); + } + } + + TYPED_TEST(fft_oracle_test, DcLandsAsNInSlotZero) { + for (const std::size_t n : closed_form_sizes()) { + std::vector x(n, 1.0); + std::vector expected(n, 0.0); + expected[0] = static_cast(n); + check_forward(x, expected, "dc"); + } + } + + TYPED_TEST(fft_oracle_test, NyquistAlternationLandsAsNInSlotOne) { + for (const std::size_t n : closed_form_sizes()) { + std::vector x(n); + for (std::size_t j = 0; j < n; ++j) { + x[j] = (j % 2 == 0) ? 1.0 : -1.0; + } + std::vector expected(n, 0.0); + expected[1] = static_cast(n); + check_forward(x, expected, "nyquist"); + } + } + + TYPED_TEST(fft_oracle_test, OnBinCosineIsPlusHalfNInTheRealSlot) { + for (const std::size_t n : closed_form_sizes()) { + const std::size_t k = tone_bin(n); + const auto x = tap::dsp::test::tone(n, static_cast(k), 1.0, 0.0); + std::vector expected(n, 0.0); + expected[2 * k] = static_cast(n) / 2.0; + check_forward(x, expected, "cosine"); + } + } + + // The sign convention: a sine at bin k lands at +N/2 in the imaginary slot + // (it would be -N/2 in the engineering convention exp(-2*pi*i/N)). + TYPED_TEST(fft_oracle_test, OnBinSineIsPlusHalfNInTheImaginarySlot) { + for (const std::size_t n : closed_form_sizes()) { + const std::size_t k = tone_bin(n); + // sin(w j) = cos(w j - pi/2) + const auto x = tap::dsp::test::tone(n, static_cast(k), 1.0, -std::numbers::pi / 2.0); + std::vector expected(n, 0.0); + expected[2 * k + 1] = static_cast(n) / 2.0; + check_forward(x, expected, "sine"); + } + } + + TYPED_TEST(fft_oracle_test, TwoToneSuperposesLinearly) { + for (const std::size_t n : closed_form_sizes()) { + if (n < 8) { + continue; // needs two distinct complex bins + } + const std::size_t k1 = tone_bin(n); + const std::size_t k2 = k1 + 1; + const double a1 = 0.75; + const double a2 = -0.25; + const auto c = tap::dsp::test::tone(n, static_cast(k1), a1, 0.0); + const auto s = tap::dsp::test::tone(n, static_cast(k2), a2, -std::numbers::pi / 2.0); + std::vector x(n); + for (std::size_t j = 0; j < n; ++j) { + x[j] = c[j] + s[j]; + } + std::vector expected(n, 0.0); + expected[2 * k1] = a1 * static_cast(n) / 2.0; + expected[2 * k2 + 1] = a2 * static_cast(n) / 2.0; + check_forward(x, expected, "two-tone"); + } + } + + // ======================================================================== + // Closed forms, inverse (unnormalized: gain N/2). + // ======================================================================== + + TYPED_TEST(fft_oracle_test, InverseOfFlatSpectrumIsHalfNImpulse) { + for (const std::size_t n : closed_form_sizes()) { + std::vector a(n, 0.0); + a[0] = 1.0; + a[1] = 1.0; + for (std::size_t k = 1; k < n / 2; ++k) { + a[2 * k] = 1.0; + } + std::vector expected(n, 0.0); + expected[0] = static_cast(n) / 2.0; + check_inverse(a, expected, "inverse impulse"); + } + } + + TYPED_TEST(fft_oracle_test, InverseOfDcSpectrumIsHalfNConstant) { + for (const std::size_t n : closed_form_sizes()) { + std::vector a(n, 0.0); + a[0] = 1.0; + std::vector expected(n, 0.5); + check_inverse(a, expected, "inverse dc"); + } + } + + TYPED_TEST(fft_oracle_test, InverseOfNyquistSpectrumIsHalfNAlternation) { + for (const std::size_t n : closed_form_sizes()) { + std::vector a(n, 0.0); + a[1] = 1.0; + std::vector expected(n); + for (std::size_t k = 0; k < n; ++k) { + expected[k] = (k % 2 == 0) ? 0.5 : -0.5; + } + check_inverse(a, expected, "inverse nyquist"); + } + } + + TYPED_TEST(fft_oracle_test, InverseOfOnBinSpectrumIsHalfNTone) { + for (const std::size_t n : closed_form_sizes()) { + const std::size_t k = tone_bin(n); + std::vector a(n, 0.0); + a[2 * k] = 1.0; // cosine part + a[2 * k + 1] = 1.0; // sine part, +i convention + // x[j] = cos(w j) + sin(w j), unnormalized inverse gain is N/2 on + // the round trip, i.e. a unit bin comes back with amplitude 1. + const auto c = tap::dsp::test::tone(n, static_cast(k), 1.0, 0.0); + const auto s = tap::dsp::test::tone(n, static_cast(k), 1.0, -std::numbers::pi / 2.0); + std::vector expected(n); + for (std::size_t j = 0; j < n; ++j) { + expected[j] = c[j] + s[j]; + } + check_inverse(a, expected, "inverse tone"); + } + } + + // ======================================================================== + // The compensated DFT, N <= 256. + // ======================================================================== + + // The oracle checks itself first: its twiddles agree with libm to double + // rounding, and its own forward followed by its own inverse is the + // identity to ~1e-15 relative, so a disagreement below is the engine's. + TEST(fft_oracle_self_check, TwiddlesMatchLibmToDoubleRounding) { + for (const std::size_t n : dft_sizes()) { + const compensated_dft oracle(n); + for (std::size_t m = 0; m < n; ++m) { + const double theta = 2.0 * std::numbers::pi * static_cast(m) / static_cast(n); + const auto& w = oracle.twiddle_at(m); + EXPECT_NEAR(dd_round(w.cos), std::cos(theta), 4.0 * std::numeric_limits::epsilon()) + << "N=" << n << " m=" << m; + EXPECT_NEAR(dd_round(w.sin), std::sin(theta), 4.0 * std::numeric_limits::epsilon()) + << "N=" << n << " m=" << m; + // The low word is a genuine correction, not noise: it is at + // most half an ulp of the high word. + EXPECT_LE(std::fabs(w.cos.lo), std::fabs(w.cos.hi) * std::numeric_limits::epsilon() + 1e-300); + } + } + } + + TEST(fft_oracle_self_check, RoundTripIsIdentityToDoubleRounding) { + for (const std::size_t n : dft_sizes()) { + const compensated_dft oracle(n); + const std::vector x = tap::dsp::test::random_signal(n, 0xC0FFEEu); + const std::vector a = oracle.forward(x); + const std::vector back = oracle.inverse_unnormalized(a); + const double gain = 2.0 / static_cast(n); + for (std::size_t i = 0; i < n; ++i) { + EXPECT_NEAR(back[i] * gain, x[i], 8.0 * std::numeric_limits::epsilon()) + << "N=" << n << " " << i; + } + } + } + + template + void expect_forward_matches_dft(std::size_t n, const std::vector& x, const char* what) { + const compensated_dft oracle(n); + const std::vector expected = oracle.forward(x); + check_forward(x, expected, what); + } + + template + void expect_inverse_matches_dft(std::size_t n, const std::vector& a, const char* what) { + const compensated_dft oracle(n); + const std::vector expected = oracle.inverse_unnormalized(a); + check_inverse(a, expected, what); + } + + TYPED_TEST(fft_oracle_test, ForwardMatchesCompensatedDftOnBroadband) { + for (const std::size_t n : dft_sizes()) { + // Drawn in the profile first so the oracle sees exactly the bits + // the engine sees; from_doubles is then the identity. + const auto x = to_doubles(tap::dsp::test::random_signal(n, 0x2545F491u)); + expect_forward_matches_dft(n, x, "dft forward broadband"); + } + } + + TYPED_TEST(fft_oracle_test, ForwardMatchesCompensatedDftOnOffBinTone) { + for (const std::size_t n : dft_sizes()) { + // Off-bin (k + 0.37) so every bin carries leakage and none is + // numerically empty: the opposite regime from the closed forms. + const auto x = + to_doubles(tap::dsp::test::tone(n, static_cast(tone_bin(n)) + 0.37, 0.8, 1.1)); + expect_forward_matches_dft(n, x, "dft forward off-bin tone"); + } + } + + TYPED_TEST(fft_oracle_test, InverseMatchesCompensatedDftOnBroadbandSpectrum) { + for (const std::size_t n : dft_sizes()) { + // An arbitrary packed spectrum, not one produced by a forward, so + // the inverse is judged on its own definition. + const auto a = to_doubles(tap::dsp::test::random_signal(n, 0x1D872B41u)); + expect_inverse_matches_dft(n, a, "dft inverse broadband"); + } + } + + TYPED_TEST(fft_oracle_test, InverseMatchesCompensatedDftOnToneSpectrum) { + for (const std::size_t n : dft_sizes()) { + const compensated_dft oracle(n); + const auto x = + to_doubles(tap::dsp::test::tone(n, static_cast(tone_bin(n)) + 0.37, 0.8, 1.1)); + // The spectrum is rounded to the profile before both sides see it. + const auto a = to_doubles(from_doubles(oracle.forward(x))); + expect_inverse_matches_dft(n, a, "dft inverse tone spectrum"); + } + } + +} // namespace diff --git a/tests/test_fft_parity_ooura.cpp b/tests/test_fft_parity_ooura.cpp new file mode 100644 index 0000000..c42ab86 --- /dev/null +++ b/tests/test_fft_parity_ooura.cpp @@ -0,0 +1,398 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// THE BIT-IDENTITY GATE FOR THE C++20 PORT OF OOURA'S rdft (Stage 2a of +// docs/audit-fft-and-code-smells.md, Part 3; the policy behind it is Part 4's +// "fp-contraction policy" section and Part 6 items N3/N4). +// +// Two sides, each re-pointable in ONE place below: +// +// ooura_ref the raw rdft / rdft_f of Takuya Ooura's fftsg.c, +// from the REFERENCE copy of the C that this target's +// CMake block compiles (tests/CMakeLists.txt). +// engine_under_test whatever the library routes the real FFT to. Today +// that is basic_real_fft, i.e. the same C, so this +// suite is Ooura-vs-Ooura and trivially green. That +// is the point: the port agent re-points this alias +// at detail::split_radix_rdft, watches the suite go +// red, and works it green statement by statement. +// +// The comparison is memcmp over the raw output bytes: not EXPECT_EQ (which +// calls +0.0 and -0.0 equal and any NaN unequal to itself), not a tolerance. +// Forward and inverse are both taken in place and unnormalized, exactly as the +// C exposes them, for every power of two from 4 to 65536 plus one run at 2^20 +// (the largest size any consumer uses; item N20), on five materials chosen to +// reach different arithmetic: broadband xorshift noise (every butterfly busy), +// an on-bin tone (one live bin, the rest numerically empty), an impulse (flat +// spectrum, no cancellation), DC (maximal cancellation in every non-zero bin) +// and a full-scale alternating +1/-1 (all energy at Nyquist). +// +// Why this target owns its compiler flags. Bit identity between a C function +// and a C++ transliteration of it holds only if both are compiled to the same +// sequence of IEEE operations. Measured for Part 4: gcc -std=c17 does not +// contract a*b+c into an FMA; gcc -std=gnu17 does; g++ contracts in both c++20 +// and gnu++20; clang contracts in every mode, statement-scoped. CMake sets none +// of this, so on any ISA with FMA (Apple arm64, Cortex-M55, x86 with -march) +// the two sides fuse differently by default and the identity is false for +// reasons that have nothing to do with the port. The gate therefore compiles +// BOTH the reference C and this TU with -ffp-contract=off (MSVC: /fp:precise +// does not contract, so nothing is passed there). Same binary, same libm: +// cross-platform identity is not claimed and does not hold today either, +// because libm's cos/sin differ in the last bit between glibc, newlib, UCRT +// and Apple (Part 4, "The float twiddles"). +// +// The same source builds a SECOND, informational target at default flags +// (TAP_DSP_PARITY_INFORMATIONAL). It prints the measured max-ulp deviation +// per N and never fails. It is not a pass/fail gate with an assumed bound +// because a different fusion choice per butterfly stage accumulates over +// log2 N stages and "1 ulp" would be a guess, not a measurement (item N3); +// its job is to put the number on the record for each platform so the policy +// fft.h states about exporting -ffp-contract=off can be decided on data. +// +// TAP_DSP_PARITY_MAX_N caps the sizes run, so the emulated QEMU legs can take +// the suite at N <= 4096 (Part 9, Part 10) where 2^20 does not fit in RAM. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "support/signals.h" +#include "tap/dsp/fft.h" + +#ifndef TAP_DSP_PARITY_MAX_N +#define TAP_DSP_PARITY_MAX_N (1 << 20) +#endif + +namespace { + + // ------------------------------------------------------------------------ + // The reference side. Re-point here, nowhere else. + // + // Raw Ooura, called exactly as fftsg.c documents: ip[0] = 0 requests table + // initialization on the first call, and the workspace geometry is the one + // readme.txt prescribes (ip: 2 + sqrt(n/2), w: n/2). rdft / rdft_f come + // from the reference C library this target links (see the CMake block), + // not from the library's tap_dsp_fft, so the reference stays the C even + // after the library stops linking it (Stage 2c). + // ------------------------------------------------------------------------ + inline void ooura_rdft_ref(int n, int isgn, double* a, int* ip, double* w) { + rdft(n, isgn, a, ip, w); + } + inline void ooura_rdft_ref(int n, int isgn, float* a, int* ip, float* w) { + rdft_f(n, isgn, a, ip, w); + } + + template + class ooura_ref { + public: + explicit ooura_ref(std::size_t n) + : m_n(static_cast(n)) + , m_ip(2 + static_cast(std::sqrt(static_cast(n) / 2.0)) + 1, 0) + , m_w(n / 2, Sample(0)) { + m_ip[0] = 0; + } + void forward_inplace(Sample* a) { ooura_rdft_ref(m_n, 1, a, m_ip.data(), m_w.data()); } + void inverse_inplace(Sample* a) { ooura_rdft_ref(m_n, -1, a, m_ip.data(), m_w.data()); } + + private: + int m_n; + std::vector m_ip; + std::vector m_w; + }; + + // ------------------------------------------------------------------------ + // The side under test. Re-point here, nowhere else. + // + // Stage 2a (the port lands beside the C, nothing routed): re-point at + // tap::dsp::detail::split_radix_rdft + // Stage 2b (routing flipped): back to basic_real_fft, which then IS the + // port, and the suite guards the flip. + // Both expose the constructor-from-size / forward_inplace / inverse_inplace + // surface Part 4 fixes, so nothing else in this file changes. + // ------------------------------------------------------------------------ + template + using engine_under_test = tap::dsp::basic_real_fft; + + // ------------------------------------------------------------------------ + // Materials. + // ------------------------------------------------------------------------ + enum class material { broadband, on_bin_tone, impulse, dc, alternating }; + + constexpr material k_materials[] = {material::broadband, material::on_bin_tone, material::impulse, material::dc, + material::alternating}; + + const char* name_of(material m) { + switch (m) { + case material::broadband: + return "broadband"; + case material::on_bin_tone: + return "on-bin tone"; + case material::impulse: + return "impulse"; + case material::dc: + return "dc"; + case material::alternating: + return "alternating +-1"; + } + return "?"; + } + + template + std::vector make_material(material m, std::size_t n) { + switch (m) { + case material::broadband: + return tap::dsp::test::random_signal(n, 0x9E3779B9u); + case material::on_bin_tone: + // Bin n/8 (bin 1 below n = 8): far from DC and Nyquist so the + // butterflies that combine it see a real twiddle, not +-1 or +-i. + return tap::dsp::test::tone(n, static_cast(std::max(1, n / 8)), 0.5, 0.3); + case material::impulse: { + std::vector x(n, Sample(0)); + x[0] = Sample(1); + return x; + } + case material::dc: + return std::vector(n, Sample(1)); + case material::alternating: { + std::vector x(n); + for (std::size_t i = 0; i < n; ++i) { + x[i] = (i % 2 == 0) ? Sample(1) : Sample(-1); + } + return x; + } + } + return {}; + } + + // ------------------------------------------------------------------------ + // Measurement. ulp distance over the ordered-integer image of the bits + // (sign-magnitude folded to two's complement), so +0/-0 count as 1 apart, + // adjacent finite values as 1, and a NaN as "far" rather than as an + // exception. Used only to make a mismatch report readable and to feed the + // informational target; the gate itself is memcmp. + // ------------------------------------------------------------------------ + template + std::uint64_t ulp_distance(Sample a, Sample b) { + using bits_t = std::conditional_t; + static_assert(sizeof(bits_t) == sizeof(Sample)); + bits_t ia = 0; + bits_t ib = 0; + std::memcpy(&ia, &a, sizeof(Sample)); + std::memcpy(&ib, &b, sizeof(Sample)); + if (ia < 0) { + ia = std::numeric_limits::min() - ia; + } + if (ib < 0) { + ib = std::numeric_limits::min() - ib; + } + const std::uint64_t ua = static_cast(ia); + const std::uint64_t ub = static_cast(ib); + return ua > ub ? ua - ub : ub - ua; + } + + struct comparison { + bool identical = true; + std::size_t first_diff = 0; + std::uint64_t max_ulp = 0; + }; + + template + comparison compare(const std::vector& got, const std::vector& ref) { + comparison c; + c.identical = std::memcmp(got.data(), ref.data(), got.size() * sizeof(Sample)) == 0; + if (c.identical) { + return c; + } + bool seen = false; + for (std::size_t i = 0; i < got.size(); ++i) { + const std::uint64_t d = ulp_distance(got[i], ref[i]); + if (d != 0 && !seen) { + c.first_diff = i; + seen = true; + } + c.max_ulp = std::max(c.max_ulp, d); + } + return c; + } + + // Runs one (N, material) pair in both directions. The inverse input is the + // reference forward's output, so both sides receive identical bits and the + // inverse is judged on its own, not through the forward. + template + struct pair_result { + comparison forward; + comparison inverse; + }; + + template + pair_result run_pair(std::size_t n, material m) { + const std::vector x = make_material(m, n); + + ooura_ref ref(n); + engine_under_test dut(n); + std::vector ref_spec = x; + std::vector dut_spec = x; + ref.forward_inplace(ref_spec.data()); + dut.forward_inplace(dut_spec.data()); + + std::vector ref_time = ref_spec; + std::vector dut_time = ref_spec; + ref.inverse_inplace(ref_time.data()); + dut.inverse_inplace(dut_time.data()); + + return {compare(dut_spec, ref_spec), compare(dut_time, ref_time)}; + } + + template + const char* precision_name() { + return sizeof(Sample) == 4 ? "float" : "double"; + } + + std::vector sizes_up_to_65536() { + std::vector sizes; + for (std::size_t n = 4; n <= 65536 && n <= static_cast(TAP_DSP_PARITY_MAX_N); n *= 2) { + sizes.push_back(n); + } + return sizes; + } + + constexpr std::size_t k_large_n = std::size_t{1} << 20; + +#if !defined(TAP_DSP_PARITY_INFORMATIONAL) + + // ======================================================================== + // THE GATE: -ffp-contract=off on both sides, memcmp identity. + // ======================================================================== + + template + void expect_identical(std::size_t n, material m, bool forward) { + const pair_result r = run_pair(n, m); + const comparison& c = forward ? r.forward : r.inverse; + ASSERT_TRUE(c.identical) << (forward ? "forward" : "inverse") << " " << precision_name() << " N=" << n + << " material=" << name_of(m) << ": first difference at index " << c.first_diff + << ", max " << c.max_ulp << " ulp. The engine under test is not the same sequence " + << "of IEEE operations as Ooura's rdft; see the statement-fidelity rule in Part 4."; + } + + template + void expect_identical_all_sizes(bool forward) { + for (const std::size_t n : sizes_up_to_65536()) { + for (const material m : k_materials) { + expect_identical(n, m, forward); + if (::testing::Test::HasFatalFailure()) { + return; + } + } + } + } + + TEST(fft_parity_ooura, ForwardIsBitIdenticalToOouraDouble) { + expect_identical_all_sizes(true); + } + + TEST(fft_parity_ooura, ForwardIsBitIdenticalToOouraFloat) { + expect_identical_all_sizes(true); + } + + TEST(fft_parity_ooura, InverseIsBitIdenticalToOouraDouble) { + expect_identical_all_sizes(false); + } + + TEST(fft_parity_ooura, InverseIsBitIdenticalToOouraFloat) { + expect_identical_all_sizes(false); + } + + // One run at 2^20 — the largest geometry any consumer uses (AmbiTap's + // long-partition convolution) — kept as its own test so its runtime shows + // separately and it can be excluded on emulated targets by name or by + // TAP_DSP_PARITY_MAX_N. + template + void expect_identical_large() { + if (k_large_n > static_cast(TAP_DSP_PARITY_MAX_N)) { + GTEST_SKIP() << "2^20 exceeds TAP_DSP_PARITY_MAX_N=" << TAP_DSP_PARITY_MAX_N; + } + for (const material m : k_materials) { + expect_identical(k_large_n, m, true); + expect_identical(k_large_n, m, false); + if (::testing::Test::HasFatalFailure()) { + return; + } + } + } + + TEST(fft_parity_ooura, LargeTransformIsBitIdenticalToOouraDouble) { + expect_identical_large(); + } + + TEST(fft_parity_ooura, LargeTransformIsBitIdenticalToOouraFloat) { + expect_identical_large(); + } + +#else // TAP_DSP_PARITY_INFORMATIONAL + + // ======================================================================== + // INFORMATIONAL: default flags on both sides, measured max ulp per N, + // never a failure. Printed to stdout so it is in every CI log, and + // recorded as gtest properties so it is in the JUnit XML too. + // ======================================================================== + + template + void report_max_ulp() { + std::vector sizes = sizes_up_to_65536(); + if (k_large_n <= static_cast(TAP_DSP_PARITY_MAX_N)) { + sizes.push_back(k_large_n); + } + std::uint64_t overall = 0; + std::printf("[ ulp ] %s: engine under test vs Ooura, default compiler flags (informational)\n", + precision_name()); + std::printf("[ ulp ] %10s %14s %14s %s\n", "N", "forward", "inverse", "worst material (fwd / inv)"); + for (const std::size_t n : sizes) { + std::uint64_t fwd = 0; + std::uint64_t inv = 0; + const char* fwd_worst = name_of(k_materials[0]); + const char* inv_worst = name_of(k_materials[0]); + for (const material m : k_materials) { + const pair_result r = run_pair(n, m); + if (r.forward.max_ulp > fwd) { + fwd = r.forward.max_ulp; + fwd_worst = name_of(m); + } + if (r.inverse.max_ulp > inv) { + inv = r.inverse.max_ulp; + inv_worst = name_of(m); + } + } + overall = std::max({overall, fwd, inv}); + std::printf("[ ulp ] %10zu %14llu %14llu %s / %s\n", n, static_cast(fwd), + static_cast(inv), fwd_worst, inv_worst); + ::testing::Test::RecordProperty(std::string(precision_name()) + "_N" + std::to_string(n) + "_fwd", + std::to_string(fwd)); + ::testing::Test::RecordProperty(std::string(precision_name()) + "_N" + std::to_string(n) + "_inv", + std::to_string(inv)); + } + std::printf("[ ulp ] %s: max over all N and materials = %llu ulp\n", precision_name(), + static_cast(overall)); + std::fflush(stdout); + SUCCEED() << "informational only: " << overall << " ulp max"; + } + + TEST(fft_parity_ooura_default_flags, ReportsMaxUlpVersusOouraDouble) { + report_max_ulp(); + } + + TEST(fft_parity_ooura_default_flags, ReportsMaxUlpVersusOouraFloat) { + report_max_ulp(); + } + +#endif // TAP_DSP_PARITY_INFORMATIONAL + +} // namespace diff --git a/tests/test_fft_rt.cpp b/tests/test_fft_rt.cpp new file mode 100644 index 0000000..018c769 --- /dev/null +++ b/tests/test_fft_rt.cpp @@ -0,0 +1,386 @@ +// SPDX-License-Identifier: MIT +// Copyright 2026 Timothy Place and the DspTap contributors. +// +// THE REAL-TIME GUARD for tap::dsp::basic_real_fft (Part 9 of +// docs/audit-fft-and-code-smells.md names this file). +// +// fft.h promises that the transforms are noexcept and allocation-free once +// the object is constructed, so they are safe on an audio thread. This file +// makes both halves of that promise checkable: +// +// - noexcept is a static_assert on the four transform entry points, on +// both instantiations (the pattern test_nn.cpp uses); +// - "allocation-free" is proved by REPLACING THE GLOBAL operator new / +// operator delete in this translation unit with counting versions and +// asserting the count does not move across a transform. The replacement +// is program-wide (that is how replaceable allocation functions work), +// so every window below is kept tight: nothing between the two reads of +// the counter but the call under test, and every buffer allocated before +// the first read. +// +// The FIRST call after construction is covered as well as a steady-state one: +// today Ooura builds its trig and bit-reversal tables lazily on the first +// transform (Part 1, item F6), and that initialization must be allocation-free +// too, because the first transform a consumer runs is very often on the audio +// thread already. The port moves table construction into the constructor +// (Part 4); this test is indifferent to where it happens, only to what it +// allocates. +// +// NOT covered, deliberately: the float-I/O convenience overloads on the +// double engine (forward(const float*, float*) / inverse(const float*, +// float*)). fft.h documents them as a setup-time path that allocates a +// staging buffer per call (Part 1, item F7) and Decision D5 deprecates them +// for one consumer cycle; guarding them would pin a behaviour the plan is +// removing. They are also not noexcept, consistently with that. +// +// Copy and copy-assignment are pinned to produce bit-identical output to the +// source in both directions and both precisions. test_fft_backend.cpp's +// CopiesAgreeWithTheirSource covers the float FORWARD at the certified +// geometries while walking the heap (its purpose is vDSP's alignment +// dispatch); this one is the general value-semantics check that the engine's +// tables travel with the object, including through operator=. +// +// Stage 4 extension points, per Part 9: the std::span overloads equal the +// pointer overloads, and Engine::is_shareable is what the header says. Neither +// exists yet, so neither is asserted here. + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "support/signals.h" +#include "tap/dsp/fft.h" + +// GCC pairs every `new` expression it inlines in this TU with the replaced +// operator delete below, sees std::free on the result of an "allocation +// function" it treats as opaque, and reports -Wmismatched-new-delete once per +// inlining site (GCC bug 101480: replaced allocation functions that forward to +// malloc/free). The pairing is correct by construction here — every allocating +// form returns std::malloc's result and every deallocating form frees it — so +// the diagnostic is silenced for this file only. +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC diagnostic ignored "-Wmismatched-new-delete" +#endif + +// ---------------------------------------------------------------------------- +// Counting replacements for the replaceable global allocation functions +// ([new.delete.single] and [new.delete.array]). All eight allocating forms +// (throwing / nothrow x plain / aligned x single / array) go through one +// counter; every deallocating form is matched so nothing falls back to the +// library default. Aligned forms are served by over-allocating and stashing +// the raw pointer just below the aligned block, which needs no platform +// aligned-malloc. +// ---------------------------------------------------------------------------- +namespace { + + std::atomic allocation_count{0}; + + void* counted_malloc(std::size_t n) { + allocation_count.fetch_add(1, std::memory_order_relaxed); + return std::malloc(n == 0 ? 1 : n); + } + + void* counted_aligned_malloc(std::size_t n, std::size_t alignment) { + allocation_count.fetch_add(1, std::memory_order_relaxed); + const std::size_t slack = alignment + sizeof(void*); + void* const raw = std::malloc(n + slack); + if (raw == nullptr) { + return nullptr; + } + const std::uintptr_t base = reinterpret_cast(raw) + sizeof(void*); + const std::uintptr_t aligned = (base + alignment - 1) & ~(static_cast(alignment) - 1); + void* const p = reinterpret_cast(aligned); + std::memcpy(reinterpret_cast(aligned) - 1, &raw, sizeof(void*)); + return p; + } + + void counted_aligned_free(void* p) noexcept { + if (p == nullptr) { + return; + } + void* raw = nullptr; + std::memcpy(&raw, static_cast(p) - 1, sizeof(void*)); + std::free(raw); + } + +} // namespace + +void* operator new(std::size_t n) { + void* p = counted_malloc(n); + if (p == nullptr) { + throw std::bad_alloc(); + } + return p; +} +void* operator new[](std::size_t n) { + void* p = counted_malloc(n); + if (p == nullptr) { + throw std::bad_alloc(); + } + return p; +} +void* operator new(std::size_t n, const std::nothrow_t&) noexcept { + return counted_malloc(n); +} +void* operator new[](std::size_t n, const std::nothrow_t&) noexcept { + return counted_malloc(n); +} +void* operator new(std::size_t n, std::align_val_t al) { + void* p = counted_aligned_malloc(n, static_cast(al)); + if (p == nullptr) { + throw std::bad_alloc(); + } + return p; +} +void* operator new[](std::size_t n, std::align_val_t al) { + void* p = counted_aligned_malloc(n, static_cast(al)); + if (p == nullptr) { + throw std::bad_alloc(); + } + return p; +} +void* operator new(std::size_t n, std::align_val_t al, const std::nothrow_t&) noexcept { + return counted_aligned_malloc(n, static_cast(al)); +} +void* operator new[](std::size_t n, std::align_val_t al, const std::nothrow_t&) noexcept { + return counted_aligned_malloc(n, static_cast(al)); +} + +void operator delete(void* p) noexcept { + std::free(p); +} +void operator delete[](void* p) noexcept { + std::free(p); +} +void operator delete(void* p, std::size_t) noexcept { + std::free(p); +} +void operator delete[](void* p, std::size_t) noexcept { + std::free(p); +} +void operator delete(void* p, const std::nothrow_t&) noexcept { + std::free(p); +} +void operator delete[](void* p, const std::nothrow_t&) noexcept { + std::free(p); +} +void operator delete(void* p, std::align_val_t) noexcept { + counted_aligned_free(p); +} +void operator delete[](void* p, std::align_val_t) noexcept { + counted_aligned_free(p); +} +void operator delete(void* p, std::size_t, std::align_val_t) noexcept { + counted_aligned_free(p); +} +void operator delete[](void* p, std::size_t, std::align_val_t) noexcept { + counted_aligned_free(p); +} +void operator delete(void* p, std::align_val_t, const std::nothrow_t&) noexcept { + counted_aligned_free(p); +} +void operator delete[](void* p, std::align_val_t, const std::nothrow_t&) noexcept { + counted_aligned_free(p); +} + +namespace { + + template + using fft_t = tap::dsp::basic_real_fft; + + // ------------------------------------------------------------------------ + // noexcept, as a compile-time fact on both instantiations. + // ------------------------------------------------------------------------ + template + constexpr bool transforms_are_noexcept() { + using fft = fft_t; + static_assert(noexcept(std::declval().forward_inplace(std::declval()))); + static_assert(noexcept(std::declval().inverse_inplace(std::declval()))); + static_assert(noexcept(std::declval().forward(std::declval(), std::declval()))); + static_assert(noexcept(std::declval().inverse(std::declval(), std::declval()))); + static_assert(noexcept(std::declval().size())); + static_assert(noexcept(std::declval().num_bins())); + return true; + } + static_assert(transforms_are_noexcept()); + static_assert(transforms_are_noexcept()); + + // ------------------------------------------------------------------------ + // The allocation guard. + // ------------------------------------------------------------------------ + class allocation_guard { + public: + allocation_guard() noexcept + : m_start(allocation_count.load(std::memory_order_relaxed)) {} + std::size_t allocations_since() const noexcept { + return allocation_count.load(std::memory_order_relaxed) - m_start; + } + + private: + std::size_t m_start; + }; + + constexpr std::size_t k_guarded_sizes[] = {512, 4096}; + + template + class fft_rt_test : public ::testing::Test {}; + + using sample_types = ::testing::Types; + TYPED_TEST_SUITE(fft_rt_test, sample_types); + + TYPED_TEST(fft_rt_test, TransformsAreNoexcept) { + EXPECT_TRUE(transforms_are_noexcept()); + } + + // The guard itself must see allocations, or every test below passes for + // the wrong reason. + TEST(fft_rt_guard, CountsAVectorAllocation) { + allocation_guard guard; + { + std::vector v(64); + EXPECT_NE(v.data(), nullptr); + } + EXPECT_GE(guard.allocations_since(), 1u); + } + + template + void expect_no_allocation(std::size_t n, const char* what, Call&& call) { + fft_t fft(n); + std::vector in = tap::dsp::test::random_signal(n, 0x9E3779B9u); + std::vector out = in; + + // First call after construction (lazy table init today, F6). + { + allocation_guard guard; + call(fft, in.data(), out.data()); + const std::size_t count = guard.allocations_since(); + EXPECT_EQ(count, 0u) << what << " N=" << n << " allocated " << count + << " time(s) on the FIRST call after construction"; + } + // Steady state. + { + allocation_guard guard; + call(fft, in.data(), out.data()); + const std::size_t count = guard.allocations_since(); + EXPECT_EQ(count, 0u) << what << " N=" << n << " allocated " << count << " time(s) in steady state"; + } + } + + TYPED_TEST(fft_rt_test, ForwardInplaceAllocatesNothing) { + for (const std::size_t n : k_guarded_sizes) { + expect_no_allocation( + n, "forward_inplace", + [](fft_t& fft, TypeParam*, TypeParam* out) { fft.forward_inplace(out); }); + } + } + + TYPED_TEST(fft_rt_test, InverseInplaceAllocatesNothing) { + for (const std::size_t n : k_guarded_sizes) { + expect_no_allocation( + n, "inverse_inplace", + [](fft_t& fft, TypeParam*, TypeParam* out) { fft.inverse_inplace(out); }); + } + } + + TYPED_TEST(fft_rt_test, ForwardOutOfPlaceAllocatesNothing) { + for (const std::size_t n : k_guarded_sizes) { + expect_no_allocation( + n, "forward", [](fft_t& fft, const TypeParam* in, TypeParam* out) { fft.forward(in, out); }); + } + } + + TYPED_TEST(fft_rt_test, InverseOutOfPlaceAllocatesNothing) { + for (const std::size_t n : k_guarded_sizes) { + expect_no_allocation( + n, "inverse", [](fft_t& fft, const TypeParam* in, TypeParam* out) { fft.inverse(in, out); }); + } + } + + // ------------------------------------------------------------------------ + // Value semantics: a copy computes exactly what its source computes. + // ------------------------------------------------------------------------ + template + void expect_bit_identical(const std::vector& a, const std::vector& b, const char* what) { + ASSERT_EQ(a.size(), b.size()); + for (std::size_t i = 0; i < a.size(); ++i) { + // memcmp, not ==: the comparison is over the exact bits. + ASSERT_EQ(std::memcmp(&a[i], &b[i], sizeof(Sample)), 0) + << what << " differs from its source at index " << i; + } + } + + template + struct both_directions { + std::vector spectrum; + std::vector time; + }; + + template + both_directions run_both(fft_t& fft, const std::vector& x) { + both_directions r; + r.spectrum = x; + fft.forward_inplace(r.spectrum.data()); + r.time = r.spectrum; + fft.inverse_inplace(r.time.data()); + return r; + } + + TYPED_TEST(fft_rt_test, CopyProducesBitIdenticalOutput) { + constexpr std::size_t n = 512; + const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); + + fft_t original(n); + // Source already warmed (tables built) so the copy carries built tables. + const auto from_original = run_both(original, x); + + fft_t copy(original); + const auto from_copy = run_both(copy, x); + expect_bit_identical(from_copy.spectrum, from_original.spectrum, "copy forward"); + expect_bit_identical(from_copy.time, from_original.time, "copy inverse"); + + // And the source is unaffected by having been copied. + const auto again = run_both(original, x); + expect_bit_identical(again.spectrum, from_original.spectrum, "source forward after copy"); + expect_bit_identical(again.time, from_original.time, "source inverse after copy"); + } + + TYPED_TEST(fft_rt_test, CopyOfUnwarmedSourceProducesBitIdenticalOutput) { + constexpr std::size_t n = 512; + const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); + + // Copied BEFORE any transform: with lazy tables (F6) both objects + // build their own; with constructor-built tables both carry the same. + fft_t original(n); + fft_t copy(original); + const auto from_copy = run_both(copy, x); + const auto from_original = run_both(original, x); + expect_bit_identical(from_copy.spectrum, from_original.spectrum, "unwarmed copy forward"); + expect_bit_identical(from_copy.time, from_original.time, "unwarmed copy inverse"); + } + + TYPED_TEST(fft_rt_test, CopyAssignmentProducesBitIdenticalOutput) { + constexpr std::size_t n = 512; + const auto x = tap::dsp::test::random_signal(n, 0x2545F491u); + + fft_t original(n); + const auto from_original = run_both(original, x); + + // Assigned over an engine of a DIFFERENT geometry, so the assignment + // has to replace every table, not just refresh one of the same size. + fft_t target(4096); + target = original; + EXPECT_EQ(target.size(), n); + EXPECT_EQ(target.num_bins(), n / 2 + 1); + const auto from_target = run_both(target, x); + expect_bit_identical(from_target.spectrum, from_original.spectrum, "assigned forward"); + expect_bit_identical(from_target.time, from_original.time, "assigned inverse"); + } + +} // namespace From fdd470d4d400f89eb4cb7e8ba1bf3721ad466c90 Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Thu, 17 Sep 2026 20:05:28 +0000 Subject: [PATCH 2/4] Address the wave-1 reviews on the Stage 2a suite (host-side items) Review findings A2-A12 and B3, B4, B6-B10 on #24; the #17-coupled items (bare-metal registration through the shared helper) follow in a second commit once the helper lands. - Informational target: a second private reference library at default flags replaces the link to tap::dsp_fft (which Stage 2c removes and which carries CMSIS objects on the M55); the report is one test for both precisions and --gtest_output=xml writes its RecordProperty values to parity-ulp.xml in the build tree; the test-file comment no longer claims the table is in every CI log (ctest hides passing stdout; CI runs the label with -V as its own step, Part 13). - TAP_DSP_PARITY_MAX_N defaults to 4096 when CMAKE_CROSSCOMPILING, is static_asserted to be a power of two >= 4, and sizes_up_to_65536() is renamed sweep_sizes(). MSVC /fp:precise claim qualified (VS 2022 17.0+). objdump FMA-count recipe and the Stage 2c fftsg_float.c hazard recorded in the parity file header. - Oracle: profile carries forward_scale(n), inverse_scale(n), tolerance(n, norm2) and k_full_scale so Q15/Q31 plug in without touching the tests; the Higham comment states mu ~ few u (makewt derives half the table arithmetically), that the constant is derived with slack (measured 0.25/0.17), and why the 2-norm-per-element bound confines the DFT comparison to N <= 256; the inverse edge term goes through TwoSum; ImpulseHasFlatSpectrum renamed ImpulseIsFlatAtEverySize and the overlap with test_fft.cpp stated in the header; self-check constants derived and their measured maxima recorded. - RT guard: the header states it counts C++ operator new only (malloc and vendor-internal allocation are invisible); the GCC pragma is push/pop scoped to the replacement functions; TransformsAreNoexcept says why it exists beside the static_asserts. - signals.h: db_from_ratio dropped until a caller exists. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- tests/CMakeLists.txt | 94 ++++++++---- tests/support/signals.h | 17 +-- tests/test_fft_oracle.cpp | 255 +++++++++++++++++++------------- tests/test_fft_parity_ooura.cpp | 56 +++++-- tests/test_fft_rt.cpp | 43 ++++-- 5 files changed, 293 insertions(+), 172 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4872cb7..53e883e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -144,7 +144,7 @@ target_link_libraries(tap_dsp_tests PRIVATE # ============================================================================== # BEGIN Stage 2a additions (docs/audit-fft-and-code-smells.md, Part 9): the # independent oracle and the real-time guard join the main executable; the -# Ooura bit-identity gate gets its own target so it can own its compiler flags. +# Ooura bit-identity gate gets its own targets so it can own its compiler flags. # Keep this block self-contained and at the end of the file. # ============================================================================== @@ -153,28 +153,53 @@ target_sources(tap_dsp_tests PRIVATE test_fft_rt.cpp) # ------------------------------------------------------------------------------ -# The reference copy of the vendored C, compiled with fp-contraction OFF. Bit -# identity between the C and its C++ transliteration holds only if both are -# compiled to the same sequence of IEEE operations, and CMake gives neither -# side an fp-contract setting: gcc contracts C in gnu17 mode, g++ contracts -# C++ in every mode, clang contracts everywhere (Part 4, "fp-contraction -# policy"). MSVC's default /fp:precise does not contract, so nothing is passed -# there. This library is private to the tests and is not the tap_dsp_fft the -# library links; the parity executable links THIS and not tap::dsp, so its -# rdft/rdft_f resolve here and no backend define (vDSP on Apple) reaches it. +# Two private reference builds of the vendored C. Neither is the tap_dsp_fft +# the library links: that target disappears at Stage 2c (it survives only under +# TAP_DSP_FFT_CMSIS, where it also carries CMSIS objects), and the gate's whole +# point is to compare against the C independently of how the library builds it. +# The parity executables link one of these and NOT tap::dsp, so their +# rdft/rdft_f resolve here and no backend define (vDSP on Apple) reaches them. +# +# Stage 2c note: fftsg.c moves to tests/reference/ooura/ (D6); this gate needs +# fftsg_float.c to move with it, or its float side becomes port-vs-port. +# +# ..._nocontract fp-contraction OFF. Bit identity between the C and its C++ +# transliteration holds only if both compile to the same +# sequence of IEEE operations, and CMake gives neither side an +# fp-contract setting: gcc contracts C in gnu17 mode, g++ +# contracts C++ in every mode, clang contracts everywhere +# (Part 4, "fp-contraction policy"). MSVC: nothing is passed; +# its default /fp:precise has not contracted since VS 2022 +# 17.0 made /fp:contract opt-in (earlier x64 /arch:AVX2 and +# ARM64 builds did contract under /fp:precise; windows-latest +# is well past 17.0). +# ..._default the same two files at default flags, for the informational +# target, so it measures "port at default flags vs C at +# default flags" and nothing else. # ------------------------------------------------------------------------------ set(_tap_dsp_nocontract_c "$<$:-ffp-contract=off>") set(_tap_dsp_nocontract_cxx "$<$:-ffp-contract=off>") - -add_library(tap_dsp_fft_reference_nocontract STATIC +set(_tap_dsp_reference_sources ${PROJECT_SOURCE_DIR}/third_party/ooura/fftsg.c ${PROJECT_SOURCE_DIR}/third_party/ooura/fftsg_float.c) + +add_library(tap_dsp_fft_reference_nocontract STATIC ${_tap_dsp_reference_sources}) target_compile_options(tap_dsp_fft_reference_nocontract PRIVATE ${_tap_dsp_nocontract_c}) -# Largest transform the parity suite runs. The default covers every consumer -# geometry (2^20); emulated targets set 4096 (Part 10) to fit in RAM. -set(TAP_DSP_PARITY_MAX_N 1048576 CACHE STRING - "Largest FFT size the Ooura parity suite runs (power of two; 4096 on emulated targets)") +add_library(tap_dsp_fft_reference_default STATIC ${_tap_dsp_reference_sources}) + +# Largest transform the parity suite runs. Hosted, the default covers every +# consumer geometry (2^20: five 8 MB buffers plus tables). Cross-compiled, the +# emulated Cortex-M legs (Part 10) get 4096, the largest size that fits the +# MPS2/MPS3 data regions; a toolchain file may still set its own cache value, +# which this default does not override. +if(CMAKE_CROSSCOMPILING) + set(_tap_dsp_parity_max_n_default 4096) +else() + set(_tap_dsp_parity_max_n_default 1048576) +endif() +set(TAP_DSP_PARITY_MAX_N ${_tap_dsp_parity_max_n_default} CACHE STRING + "Largest FFT size the Ooura parity suite runs (power of two >= 4; 4096 on emulated targets)") # THE GATE: -ffp-contract=off on both sides, memcmp bit identity. add_executable(tap_dsp_fft_parity test_fft_parity_ooura.cpp) @@ -184,16 +209,18 @@ target_compile_options(tap_dsp_fft_parity PRIVATE ${_tap_dsp_nocontract_cxx}) target_compile_definitions(tap_dsp_fft_parity PRIVATE TAP_DSP_PARITY_MAX_N=${TAP_DSP_PARITY_MAX_N}) target_link_libraries(tap_dsp_fft_parity PRIVATE tap_dsp_fft_reference_nocontract - tap_dsp_warnings - GTest::gtest_main) -gtest_discover_tests(tap_dsp_fft_parity DISCOVERY_TIMEOUT 120 PROPERTIES TIMEOUT 900 LABELS parity) + tap_dsp_warnings) -# INFORMATIONAL: the same source at default flags on both sides (the C here is -# the library's own tap_dsp_fft build), printing the measured max-ulp deviation -# per N and never failing. It is not a pass/fail gate with an assumed bound -# because a different fusion choice per butterfly stage accumulates over -# log2 N stages; the number goes on the record per platform instead (Part 6, -# item N3). +# INFORMATIONAL: the same source at default flags on both sides, measuring the +# max-ulp deviation per N and never failing. It is not a pass/fail gate with an +# assumed bound because a different fusion choice per butterfly stage +# accumulates over log2 N stages; the number goes on the record per platform +# instead (Part 6, item N3). Two channels carry it, because ctest hides the +# stdout of a passing test under --output-on-failure: (1) CI runs this label +# with -V as its own step (Part 13), (2) --gtest_output=xml writes the +# RecordProperty values to parity-ulp.xml in the build tree for upload. The +# report is ONE test covering both precisions so a single invocation produces +# the complete XML (each ctest entry is a separate process writing that file). add_executable(tap_dsp_fft_parity_default_flags test_fft_parity_ooura.cpp) target_include_directories(tap_dsp_fft_parity_default_flags PRIVATE ${PROJECT_SOURCE_DIR}/include) target_compile_features(tap_dsp_fft_parity_default_flags PRIVATE cxx_std_20) @@ -201,10 +228,19 @@ target_compile_definitions(tap_dsp_fft_parity_default_flags PRIVATE TAP_DSP_PARITY_INFORMATIONAL TAP_DSP_PARITY_MAX_N=${TAP_DSP_PARITY_MAX_N}) target_link_libraries(tap_dsp_fft_parity_default_flags PRIVATE - tap::dsp_fft - tap_dsp_warnings - GTest::gtest_main) -gtest_discover_tests(tap_dsp_fft_parity_default_flags DISCOVERY_TIMEOUT 120 PROPERTIES TIMEOUT 900 LABELS parity) + tap_dsp_fft_reference_default + tap_dsp_warnings) + +include(GoogleTest) +target_link_libraries(tap_dsp_fft_parity PRIVATE GTest::gtest_main) +gtest_discover_tests(tap_dsp_fft_parity + DISCOVERY_TIMEOUT 120 + PROPERTIES TIMEOUT 900 LABELS parity) +target_link_libraries(tap_dsp_fft_parity_default_flags PRIVATE GTest::gtest_main) +gtest_discover_tests(tap_dsp_fft_parity_default_flags + EXTRA_ARGS --gtest_output=xml:${CMAKE_BINARY_DIR}/parity-ulp.xml + DISCOVERY_TIMEOUT 120 + PROPERTIES TIMEOUT 900 LABELS parity) # ============================================================================== # END Stage 2a additions diff --git a/tests/support/signals.h b/tests/support/signals.h index 6f6963d..0a6474d 100644 --- a/tests/support/signals.h +++ b/tests/support/signals.h @@ -3,12 +3,12 @@ // SPDX-License-Identifier: MIT // Copyright 2026 Timothy Place and the DspTap contributors. // -// One xorshift32, one random_signal, one tone synthesizer and one dB helper, -// so a new test file does not grow its own copy (three already exist in -// test_fft.cpp, test_fft_backend.cpp and test_nn.cpp; migrating them is the -// Stage 6 hygiene item in docs/audit-fft-and-code-smells.md, not this file's -// job). Everything here is deterministic — fixed seeds, no wall clock, no -// filesystem — so it can run unchanged on the bare-metal QEMU legs. +// One xorshift32, one random_signal and one tone synthesizer, so a new test +// file does not grow its own copy (three already exist in test_fft.cpp, +// test_fft_backend.cpp and test_nn.cpp; migrating them, and adding the dB +// helper Part 9 lists once a caller exists, is the Stage 6 hygiene item in +// docs/audit-fft-and-code-smells.md, not this file's job). Everything here is deterministic — fixed seeds, no wall +// clock, no filesystem — so it can run unchanged on the bare-metal QEMU legs. #pragma once @@ -75,9 +75,4 @@ namespace tap::dsp::test { return x; } - /// 20*log10(ratio) of an amplitude ratio; -inf for zero, as std::log10 gives. - inline double db_from_ratio(double ratio) noexcept { - return 20.0 * std::log10(ratio); - } - } // namespace tap::dsp::test diff --git a/tests/test_fft_oracle.cpp b/tests/test_fft_oracle.cpp index 754da1e..e4b4fca 100644 --- a/tests/test_fft_oracle.cpp +++ b/tests/test_fft_oracle.cpp @@ -16,7 +16,7 @@ // W = exp(+2*pi*i/N) convention, the conjugate of the engineering DFT), // and a two-tone superposition — and their unnormalized inverses, whose // gain is N/2 per the contract in fft.h. Checked to a tolerance DERIVED -// from the sample type's epsilon and log2 N (see tolerance() below). +// from the sample type's epsilon and log2 N (see "Tolerance" below). // // 2. A compensated-summation DFT for N <= 256: the O(N^2) definition with // double-double accumulation (TwoSum / TwoProd, Dekker 1971 and Knuth @@ -26,11 +26,19 @@ // `double` on MSVC and on Apple arm64 (Part 6, item N19), so it would be // no oracle at all on two of the three hosted CI legs. // +// Relation to test_fft.cpp. That file's contract battery already pins the +// impulse, DC/Nyquist packing and the +i sign convention at one size each +// (N = 64 and 128) with a round tolerance; Part 9 lists the closed forms in +// both files on purpose. What this file adds is the sweep to N = 65536, the +// inverse closed forms, the derived (not round) tolerance, and the DFT +// reference; the test names are distinct so the ctest listing carries each +// promise once. +// // The suite is typed over float and double now. The fixed-point stage -// (Stage 3, Part 7) extends `profile` below with int16_t and int32_t, -// supplying the profile's documented forward scale and the conversion into -// the double domain; the tests themselves are written against that trait so -// nothing else changes. +// (Stage 3, Part 7) extends `profile` below with int16_t and int32_t; +// the tests are written against that trait (scale as a function of N, the +// profile's own tolerance, a full-scale amplitude below saturation), so the +// change is confined to the trait and the type list. #include #include @@ -47,83 +55,96 @@ namespace { + // ------------------------------------------------------------------------ + // Tolerance. + // + // Higham, Accuracy and Stability of Numerical Algorithms (2nd ed., 2002), + // Theorem 24.2: a radix-2 FFT of length N = 2^t whose twiddles carry + // relative error at most mu computes y_hat with + // + // || y_hat - y ||_2 <= t * eta / (1 - t * eta) * || y ||_2, + // eta = mu + gamma_4 (1 + mu), gamma_4 ~ 4u, + // + // where u is the unit roundoff (epsilon / 2). mu is NOT one rounding here: + // Ooura's makewt takes cos/sin from libm for a quarter of the table and + // derives the rest arithmetically (w[2] = 0.5 / cos(2 delta), the halving + // recurrences 0.5 / wk1r, ...), and the float instantiation forms + // delta * j in float, so mu is a small multiple of u. With mu ~ 2-3u, + // eta ~ 6-7u and, for t * eta << 1, the 2-norm error is below about + // 7 * u * log2(N) * ||y||_2. The largest single-element error cannot exceed + // the 2-norm, and ||y||_2 comes free with the exact answer. Casting the + // closed-form input to Sample adds at most u * ||y||_2 (Parseval), and + // Ooura is split-radix rather than the theorem's radix-2, which changes the + // constant by a factor of order one. The bound is taken as + // + // k_higham_constant * epsilon * log2(N) * ||y||_2, k_higham_constant = 4 + // + // (epsilon = 2u, so 8u: the ~7u above plus the input cast), i.e. derived + // with slack, not fitted. Measured against it on x86-64 Linux (GCC 13, -O3, + // glibc, Ooura as the engine): the largest |error| / (epsilon * log2(N) * + // ||y||_2) over every test in this file is 0.25 for double and 0.17 for + // float, 16x and 23x inside the constant. + // + // What the bound cannot see. A 2-norm used per element is loose by up to + // sqrt(N): at float, N = 65536, DC input, the tolerance is ~0.5 on a + // spectrum whose one live bin is 65536, so a single empty bin that is off + // by 0.4 would pass. That is acceptable for the closed forms (their job is + // the sweep over sizes and the sign convention) and is exactly why the DFT + // comparison is confined to N <= 256, where sqrt(N) <= 16 and the reference + // is exact to ~1e-30: there the check is tight enough to catch a single + // wrong bin. Precision pins as measured numbers live in test_fft.cpp and + // test_fft_backend.cpp; this file is an oracle for correctness. + // ------------------------------------------------------------------------ + constexpr double k_higham_constant = 4.0; + + double higham_tolerance(double epsilon, std::size_t n, double spectrum_norm2) { + const double log2n = std::log2(static_cast(n)); + return k_higham_constant * epsilon * log2n * spectrum_norm2; + } + // ------------------------------------------------------------------------ // Per-profile traits — THE EXTENSION POINT FOR THE FIXED-POINT STAGE. // - // A profile says how to get a sample into the double domain, what its - // unit roundoff is, and what scale the engine's forward / unnormalized - // inverse carry relative to the mathematical DFT (1 and 1 for the float - // profiles; the Q15/Q31 profiles document X/N under `fixed` scaling and - // will state their own inverse scale). Stage 3 adds: + // A profile says how a sample crosses into the double domain, the amplitude + // the closed forms are driven at, what scale the engine's forward and + // unnormalized inverse carry relative to the mathematical DFT as a function + // of N, and its own tolerance. For the float profiles: full scale 1.0, both + // scales 1, and the Higham bound above. Stage 3 (Part 7) adds // template <> struct profile { ... }; // template <> struct profile { ... }; - // and appends the two types to oracle_types. The tolerance for those - // profiles is a quantization-noise number from Part 7, not epsilon-based; - // tolerance() dispatches on the trait so that is a local change too. + // with k_full_scale below 1 - 2^-15 (a full-scale 1.0 input whose X/N + // expectation is exactly 1.0 saturates in Q15), forward_scale(n) = 1/n + // under `fixed` scaling (BFP reports its exponent alongside), the profile's + // inverse scale, and a tolerance built from Part 7's quantization-noise + // numbers rather than from epsilon. Then append the types to oracle_types. // ------------------------------------------------------------------------ template struct profile; template <> struct profile { - static constexpr double k_epsilon = std::numeric_limits::epsilon(); - static constexpr double k_forward_scale = 1.0; ///< engine forward = k * DFT - static constexpr double k_inverse_scale = 1.0; ///< engine inverse = k * Ooura's unnormalized inverse + static constexpr double k_epsilon = std::numeric_limits::epsilon(); + static constexpr double k_full_scale = 1.0; ///< closed-form drive amplitude static double to_double(double v) { return v; } static double from_double(double v) { return v; } + static double forward_scale(std::size_t) { return 1.0; } ///< engine forward = scale * DFT + static double inverse_scale(std::size_t) { return 1.0; } ///< engine inverse = scale * unnormalized + static double tolerance(std::size_t n, double norm2) { return higham_tolerance(k_epsilon, n, norm2); } }; template <> struct profile { - static constexpr double k_epsilon = std::numeric_limits::epsilon(); - static constexpr double k_forward_scale = 1.0; - static constexpr double k_inverse_scale = 1.0; + static constexpr double k_epsilon = std::numeric_limits::epsilon(); + static constexpr double k_full_scale = 1.0; static double to_double(float v) { return static_cast(v); } static float from_double(double v) { return static_cast(v); } + static double forward_scale(std::size_t) { return 1.0; } + static double inverse_scale(std::size_t) { return 1.0; } + static double tolerance(std::size_t n, double norm2) { return higham_tolerance(k_epsilon, n, norm2); } }; using oracle_types = ::testing::Types; - // ------------------------------------------------------------------------ - // Tolerance. - // - // Higham, Accuracy and Stability of Numerical Algorithms (2nd ed., 2002), - // Theorem 24.2: a radix-2 FFT of length N = 2^t whose twiddles carry - // relative error at most mu computes y_hat with - // - // || y_hat - y ||_2 <= t * eta / (1 - t * eta) * || y ||_2, - // eta = mu + gamma_4 (1 + mu), gamma_4 ~ 4u, - // - // where u is the unit roundoff (epsilon / 2). Ooura's tables are computed - // from libm and stored in Sample, so mu <= u, hence eta <= 5u + O(u^2) - // and, for t*eta << 1, the 2-norm error is below 5 * u * log2(N) * ||y||_2. - // The largest single-element error cannot exceed the 2-norm, and a - // measured (not remembered) value for || y ||_2 comes free with the exact - // answer. Two further terms ride on top: casting the closed-form input to - // Sample perturbs each x_j by at most u|x_j|, which by Parseval moves y by - // at most u * ||y||_2 in 2-norm; and Ooura is split-radix, not the radix-2 - // of the theorem, so its constant differs in a factor of order one. The - // bound is therefore taken as k_higham_constant * epsilon * log2(N) * - // ||y||_2 with k_higham_constant = 4 (epsilon = 2u, so 4 * epsilon = 8u: - // the 5u of the theorem, the u of the input cast, and 2u of margin for - // the split-radix constant). - // - // Measured against that bound on x86-64 Linux (GCC 13, -O3, glibc libm, - // Ooura as the engine): the largest ratio |error| / (epsilon * log2(N) * - // ||y||_2) over every test in this file was 0.25 for double and 0.17 for - // float, i.e. the engine sits 16x (double) and 23x (float) inside the - // derived bound of 4. That margin is deliberate: this is an oracle for - // correctness, not a precision ratchet; the measured-number pins live in - // test_fft.cpp and test_fft_backend.cpp. - // ------------------------------------------------------------------------ - constexpr double k_higham_constant = 4.0; - - template - double tolerance(std::size_t n, double spectrum_norm2) { - const double log2n = std::log2(static_cast(n)); - return k_higham_constant * profile::k_epsilon * log2n * spectrum_norm2; - } - /// ||y||_2 of the full complex spectrum of x, from Parseval: sqrt(N) * ||x||_2. double spectrum_norm2(const std::vector& x) { double e = 0.0; @@ -133,6 +154,13 @@ namespace { return std::sqrt(static_cast(x.size()) * e); } + std::vector scale_by(std::vector x, double factor) { + for (double& v : x) { + v *= factor; + } + return x; + } + template std::vector to_doubles(const std::vector& x) { std::vector d(x.size()); @@ -336,8 +364,9 @@ namespace { std::vector inverse_unnormalized(const std::vector& a) const { std::vector x(m_n, 0.0); for (std::size_t k = 0; k < m_n; ++k) { - const double edge = (k % 2 == 0) ? a[0] + a[1] : a[0] - a[1]; - dd acc = dd_mul_d({edge, 0.0}, 0.5); + // (R[0] +- R[N/2]) / 2: the sum error-free via TwoSum, the + // halving exact. + dd acc = dd_mul_d(two_sum(a[0], (k % 2 == 0) ? a[1] : -a[1]), 0.5); for (std::size_t j = 1; j < m_n / 2; ++j) { const dd_sincos& w = m_twiddles[(j * k) % m_n]; acc = dd_add(acc, dd_mul_d(w.cos, a[2 * j])); @@ -391,32 +420,28 @@ namespace { } // Runs the engine forward on x (given in the double domain, cast to the - // profile) and checks it against the exact spectrum `expected`. + // profile) and checks it against the exact spectrum `expected`, with the + // profile's forward scale applied. template void check_forward(const std::vector& x, const std::vector& expected, const char* what) { const std::size_t n = x.size(); tap::dsp::basic_real_fft fft(n); std::vector buf = from_doubles(x); fft.forward_inplace(buf.data()); - std::vector scaled = expected; - for (double& v : scaled) { - v *= profile::k_forward_scale; - } - expect_close(buf, scaled, tolerance(n, spectrum_norm2(x)), what, n); + const std::vector scaled = scale_by(expected, profile::forward_scale(n)); + expect_close(buf, scaled, profile::tolerance(n, spectrum_norm2(x)), what, n); } // Runs the engine's UNNORMALIZED inverse on the packed spectrum `a` and - // checks it against `expected` (which already carries the N/2 gain). + // checks it against `expected` (which already carries the N/2 gain), with + // the profile's inverse scale applied. template void check_inverse(const std::vector& a, const std::vector& expected, const char* what) { const std::size_t n = a.size(); tap::dsp::basic_real_fft fft(n); std::vector buf = from_doubles(a); fft.inverse_inplace(buf.data()); - std::vector scaled = expected; - for (double& v : scaled) { - v *= profile::k_inverse_scale; - } + const std::vector scaled = scale_by(expected, profile::inverse_scale(n)); // ||expected||_2 plays the role of ||y||_2 for the inverse direction; // the packed spectrum's 2-norm is within sqrt(2) of the true one, and // the N/2 gain means ||expected||_2 >= ||a||_2 * sqrt(N)/2 for these @@ -425,7 +450,7 @@ namespace { for (const double v : expected) { e += v * v; } - expect_close(buf, scaled, tolerance(n, std::sqrt(e)), what, n); + expect_close(buf, scaled, profile::tolerance(n, std::sqrt(e)), what, n); } template @@ -433,50 +458,55 @@ namespace { TYPED_TEST_SUITE(fft_oracle_test, oracle_types); // ======================================================================== - // Closed forms, forward. + // Closed forms, forward. Every one is driven at the profile's full-scale + // amplitude (1.0 for float and double); the exact answers scale with it. // ======================================================================== - TYPED_TEST(fft_oracle_test, ImpulseHasFlatSpectrum) { + TYPED_TEST(fft_oracle_test, ImpulseIsFlatAtEverySize) { + const double a = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { std::vector x(n, 0.0); - x[0] = 1.0; + x[0] = a; std::vector expected(n, 0.0); - expected[0] = 1.0; // DC - expected[1] = 1.0; // Nyquist + expected[0] = a; // DC + expected[1] = a; // Nyquist for (std::size_t k = 1; k < n / 2; ++k) { - expected[2 * k] = 1.0; + expected[2 * k] = a; } check_forward(x, expected, "impulse"); } } TYPED_TEST(fft_oracle_test, DcLandsAsNInSlotZero) { + const double a = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { - std::vector x(n, 1.0); + std::vector x(n, a); std::vector expected(n, 0.0); - expected[0] = static_cast(n); + expected[0] = a * static_cast(n); check_forward(x, expected, "dc"); } } TYPED_TEST(fft_oracle_test, NyquistAlternationLandsAsNInSlotOne) { + const double a = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { std::vector x(n); for (std::size_t j = 0; j < n; ++j) { - x[j] = (j % 2 == 0) ? 1.0 : -1.0; + x[j] = (j % 2 == 0) ? a : -a; } std::vector expected(n, 0.0); - expected[1] = static_cast(n); + expected[1] = a * static_cast(n); check_forward(x, expected, "nyquist"); } } TYPED_TEST(fft_oracle_test, OnBinCosineIsPlusHalfNInTheRealSlot) { + const double a = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { const std::size_t k = tone_bin(n); - const auto x = tap::dsp::test::tone(n, static_cast(k), 1.0, 0.0); + const auto x = tap::dsp::test::tone(n, static_cast(k), a, 0.0); std::vector expected(n, 0.0); - expected[2 * k] = static_cast(n) / 2.0; + expected[2 * k] = a * static_cast(n) / 2.0; check_forward(x, expected, "cosine"); } } @@ -484,25 +514,27 @@ namespace { // The sign convention: a sine at bin k lands at +N/2 in the imaginary slot // (it would be -N/2 in the engineering convention exp(-2*pi*i/N)). TYPED_TEST(fft_oracle_test, OnBinSineIsPlusHalfNInTheImaginarySlot) { + const double a = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { const std::size_t k = tone_bin(n); // sin(w j) = cos(w j - pi/2) - const auto x = tap::dsp::test::tone(n, static_cast(k), 1.0, -std::numbers::pi / 2.0); + const auto x = tap::dsp::test::tone(n, static_cast(k), a, -std::numbers::pi / 2.0); std::vector expected(n, 0.0); - expected[2 * k + 1] = static_cast(n) / 2.0; + expected[2 * k + 1] = a * static_cast(n) / 2.0; check_forward(x, expected, "sine"); } } TYPED_TEST(fft_oracle_test, TwoToneSuperposesLinearly) { + const double a = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { if (n < 8) { continue; // needs two distinct complex bins } const std::size_t k1 = tone_bin(n); const std::size_t k2 = k1 + 1; - const double a1 = 0.75; - const double a2 = -0.25; + const double a1 = 0.75 * a; + const double a2 = -0.25 * a; const auto c = tap::dsp::test::tone(n, static_cast(k1), a1, 0.0); const auto s = tap::dsp::test::tone(n, static_cast(k2), a2, -std::numbers::pi / 2.0); std::vector x(n); @@ -521,50 +553,55 @@ namespace { // ======================================================================== TYPED_TEST(fft_oracle_test, InverseOfFlatSpectrumIsHalfNImpulse) { + const double amp = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { std::vector a(n, 0.0); - a[0] = 1.0; - a[1] = 1.0; + a[0] = amp; + a[1] = amp; for (std::size_t k = 1; k < n / 2; ++k) { - a[2 * k] = 1.0; + a[2 * k] = amp; } std::vector expected(n, 0.0); - expected[0] = static_cast(n) / 2.0; + expected[0] = amp * static_cast(n) / 2.0; check_inverse(a, expected, "inverse impulse"); } } TYPED_TEST(fft_oracle_test, InverseOfDcSpectrumIsHalfNConstant) { + const double amp = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { std::vector a(n, 0.0); - a[0] = 1.0; - std::vector expected(n, 0.5); + a[0] = amp; + std::vector expected(n, 0.5 * amp); check_inverse(a, expected, "inverse dc"); } } TYPED_TEST(fft_oracle_test, InverseOfNyquistSpectrumIsHalfNAlternation) { + const double amp = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { std::vector a(n, 0.0); - a[1] = 1.0; + a[1] = amp; std::vector expected(n); for (std::size_t k = 0; k < n; ++k) { - expected[k] = (k % 2 == 0) ? 0.5 : -0.5; + expected[k] = (k % 2 == 0) ? 0.5 * amp : -0.5 * amp; } check_inverse(a, expected, "inverse nyquist"); } } TYPED_TEST(fft_oracle_test, InverseOfOnBinSpectrumIsHalfNTone) { + const double amp = profile::k_full_scale; for (const std::size_t n : closed_form_sizes()) { const std::size_t k = tone_bin(n); std::vector a(n, 0.0); - a[2 * k] = 1.0; // cosine part - a[2 * k + 1] = 1.0; // sine part, +i convention - // x[j] = cos(w j) + sin(w j), unnormalized inverse gain is N/2 on - // the round trip, i.e. a unit bin comes back with amplitude 1. - const auto c = tap::dsp::test::tone(n, static_cast(k), 1.0, 0.0); - const auto s = tap::dsp::test::tone(n, static_cast(k), 1.0, -std::numbers::pi / 2.0); + a[2 * k] = amp; // cosine part + a[2 * k + 1] = amp; // sine part, +i convention + // x[j] = amp * (cos(w j) + sin(w j)): the unnormalized inverse has + // gain N/2 on the round trip, so a bin of amp comes back with + // amplitude amp. + const auto c = tap::dsp::test::tone(n, static_cast(k), amp, 0.0); + const auto s = tap::dsp::test::tone(n, static_cast(k), amp, -std::numbers::pi / 2.0); std::vector expected(n); for (std::size_t j = 0; j < n; ++j) { expected[j] = c[j] + s[j]; @@ -579,7 +616,17 @@ namespace { // The oracle checks itself first: its twiddles agree with libm to double // rounding, and its own forward followed by its own inverse is the - // identity to ~1e-15 relative, so a disagreement below is the engine's. + // identity to double rounding, so a disagreement below is the engine's. + // The two constants are derived bounds, not fitted numbers. Twiddles: the + // comparison forms theta = 2*pi*m/n in double, and rounding theta (at most + // half an ulp of 2*pi, 4.4e-16 = 2 eps) moves cos/sin by up to that much; + // libm adds at most 1 ulp (1 eps for values in [0.5, 1]); dd_round adds + // half an ulp. Bound 3.5 eps, constant 4 eps; measured maximum 6.4e-16 = + // 2.9 eps (x86-64 Linux, glibc 2.39), i.e. the dd twiddles are as good as + // libm and the argument rounding dominates. Round trip: forward then + // inverse re-rounds through two dd_round steps plus the final 2/N + // multiply (exact, power of two), so ~1.5 eps relative on unit-scale data; + // constant 8 eps; measured maximum 1.1e-16 = 0.5 eps. TEST(fft_oracle_self_check, TwiddlesMatchLibmToDoubleRounding) { for (const std::size_t n : dft_sizes()) { const compensated_dft oracle(n); diff --git a/tests/test_fft_parity_ooura.cpp b/tests/test_fft_parity_ooura.cpp index c42ab86..1dfd9eb 100644 --- a/tests/test_fft_parity_ooura.cpp +++ b/tests/test_fft_parity_ooura.cpp @@ -41,16 +41,38 @@ // because libm's cos/sin differ in the last bit between glibc, newlib, UCRT // and Apple (Part 4, "The float twiddles"). // +// To see that the flag is load-bearing rather than take it on faith, count +// the fused instructions in the reference object on an FMA-capable target: +// +// cc -O3 -march=haswell [-ffp-contract=off] -c third_party/ooura/fftsg.c -o probe.o +// objdump -d probe.o | grep -ciE 'vfmadd|vfmsub|vfnmadd|vfnmsub' +// +// Measured: gcc 13.3 200 / clang 18.1 171 fused instructions by default, 0 +// with the flag, for both compilers. On plain x86-64 (no -march) the count is +// 0 either way because the ISA has no FMA, which is why the hosted CI legs +// cannot see a missing flag; the objdump check is the one that can. +// // The same source builds a SECOND, informational target at default flags -// (TAP_DSP_PARITY_INFORMATIONAL). It prints the measured max-ulp deviation -// per N and never fails. It is not a pass/fail gate with an assumed bound -// because a different fusion choice per butterfly stage accumulates over -// log2 N stages and "1 ulp" would be a guess, not a measurement (item N3); -// its job is to put the number on the record for each platform so the policy -// fft.h states about exporting -ffp-contract=off can be decided on data. +// (TAP_DSP_PARITY_INFORMATIONAL). It measures the max-ulp deviation per N and +// never fails. It is not a pass/fail gate with an assumed bound because a +// different fusion choice per butterfly stage accumulates over log2 N stages +// and "1 ulp" would be a guess, not a measurement (item N3); its job is to +// put the number on the record for each platform so the policy fft.h states +// about exporting -ffp-contract=off can be decided on data. Where the number +// lands: ctest hides the stdout of a passing test under --output-on-failure, +// so the table is visible only when the label is run with -V (CI does that +// as its own step, Part 13) and, as RecordProperty values, in the JUnit XML +// that the CMake block requests with --gtest_output=xml. // // TAP_DSP_PARITY_MAX_N caps the sizes run, so the emulated QEMU legs can take -// the suite at N <= 4096 (Part 9, Part 10) where 2^20 does not fit in RAM. +// the suite at N <= 4096 (Part 9, Part 10) where 2^20 does not fit in RAM; +// the CMake block defaults it to 4096 when cross-compiling. +// +// Stage 2c note. fftsg.c moves to tests/reference/ooura/ (Decision D6), and +// the 2c text deletes fftsg_float.c. This gate needs BOTH files: the float +// side compares against rdft_f, and without fftsg_float.c it degenerates to +// port-vs-port. 2c must move fftsg_float.c alongside fftsg.c, or delete the +// float half of this file and say so in the 2c PR. #include #include @@ -71,6 +93,8 @@ #ifndef TAP_DSP_PARITY_MAX_N #define TAP_DSP_PARITY_MAX_N (1 << 20) #endif +static_assert(TAP_DSP_PARITY_MAX_N >= 4 && (TAP_DSP_PARITY_MAX_N & (TAP_DSP_PARITY_MAX_N - 1)) == 0, + "TAP_DSP_PARITY_MAX_N must be a power of two >= 4, or every sweep below is vacuously green"); namespace { @@ -257,7 +281,8 @@ namespace { return sizeof(Sample) == 4 ? "float" : "double"; } - std::vector sizes_up_to_65536() { + /// The sweep: every power of two from 4 to min(65536, TAP_DSP_PARITY_MAX_N). + std::vector sweep_sizes() { std::vector sizes; for (std::size_t n = 4; n <= 65536 && n <= static_cast(TAP_DSP_PARITY_MAX_N); n *= 2) { sizes.push_back(n); @@ -285,7 +310,7 @@ namespace { template void expect_identical_all_sizes(bool forward) { - for (const std::size_t n : sizes_up_to_65536()) { + for (const std::size_t n : sweep_sizes()) { for (const material m : k_materials) { expect_identical(n, m, forward); if (::testing::Test::HasFatalFailure()) { @@ -341,13 +366,15 @@ namespace { // ======================================================================== // INFORMATIONAL: default flags on both sides, measured max ulp per N, - // never a failure. Printed to stdout so it is in every CI log, and - // recorded as gtest properties so it is in the JUnit XML too. + // never a failure. Printed to stdout (visible under `ctest -V`, NOT under + // --output-on-failure, which hides a passing test's output) and recorded + // as gtest properties for the JUnit XML the CMake block requests. One + // test for both precisions, so one process writes the whole XML. // ======================================================================== template void report_max_ulp() { - std::vector sizes = sizes_up_to_65536(); + std::vector sizes = sweep_sizes(); if (k_large_n <= static_cast(TAP_DSP_PARITY_MAX_N)) { sizes.push_back(k_large_n); } @@ -385,11 +412,8 @@ namespace { SUCCEED() << "informational only: " << overall << " ulp max"; } - TEST(fft_parity_ooura_default_flags, ReportsMaxUlpVersusOouraDouble) { + TEST(fft_parity_ooura_default_flags, ReportsMaxUlpVersusOoura) { report_max_ulp(); - } - - TEST(fft_parity_ooura_default_flags, ReportsMaxUlpVersusOouraFloat) { report_max_ulp(); } diff --git a/tests/test_fft_rt.cpp b/tests/test_fft_rt.cpp index 018c769..3258bb8 100644 --- a/tests/test_fft_rt.cpp +++ b/tests/test_fft_rt.cpp @@ -10,7 +10,7 @@ // // - noexcept is a static_assert on the four transform entry points, on // both instantiations (the pattern test_nn.cpp uses); -// - "allocation-free" is proved by REPLACING THE GLOBAL operator new / +// - "allocation-free" is checked by REPLACING THE GLOBAL operator new / // operator delete in this translation unit with counting versions and // asserting the count does not move across a transform. The replacement // is program-wide (that is how replaceable allocation functions work), @@ -18,6 +18,18 @@ // the counter but the call under test, and every buffer allocated before // the first read. // +// WHAT THE GUARD SEES AND WHAT IT DOES NOT. It counts C++ allocation +// functions only. An engine that allocates through malloc/calloc directly, or +// inside a vendor library (Apple's vDSP on the macOS leg, CMSIS on the M55), +// is invisible to it. For today's Ooura path the claim is complete: makewt +// and makect write into the caller's preallocated ip/w tables (fftsg.c +// 655-756) and the C never calls an allocator. For the backends the guard +// covers the wrapper's own code and nothing more; a malloc interposer +// (glibc's __libc_malloc, or DYLD_INTERPOSE) is the tool for the vendor +// layer and is out of scope here. The self-test CountsAVectorAllocation +// keeps the counter honest: if the replacement ever stopped being picked up, +// that test fails first. +// // The FIRST call after construction is covered as well as a steady-state one: // today Ooura builds its trig and bit-reversal tables lazily on the first // transform (Part 1, item F6), and that initialization must be allocation-free @@ -58,17 +70,6 @@ #include "support/signals.h" #include "tap/dsp/fft.h" -// GCC pairs every `new` expression it inlines in this TU with the replaced -// operator delete below, sees std::free on the result of an "allocation -// function" it treats as opaque, and reports -Wmismatched-new-delete once per -// inlining site (GCC bug 101480: replaced allocation functions that forward to -// malloc/free). The pairing is correct by construction here — every allocating -// form returns std::malloc's result and every deallocating form frees it — so -// the diagnostic is silenced for this file only. -#if defined(__GNUC__) && !defined(__clang__) -#pragma GCC diagnostic ignored "-Wmismatched-new-delete" -#endif - // ---------------------------------------------------------------------------- // Counting replacements for the replaceable global allocation functions // ([new.delete.single] and [new.delete.array]). All eight allocating forms @@ -112,6 +113,18 @@ namespace { } // namespace +// GCC pairs the `new` expressions it inlines in this TU with the replaced +// operator delete below, sees std::free on the result of an "allocation +// function" it treats as opaque, and reports -Wmismatched-new-delete once per +// inlining site (GCC bug 101480: replaced allocation functions that forward to +// malloc/free). The pairing is correct by construction here — every allocating +// form returns std::malloc's result and every deallocating form frees it — so +// the diagnostic is silenced around the replacement functions only. +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wmismatched-new-delete" +#endif + void* operator new(std::size_t n) { void* p = counted_malloc(n); if (p == nullptr) { @@ -190,6 +203,10 @@ void operator delete[](void* p, std::align_val_t, const std::nothrow_t&) noexcep counted_aligned_free(p); } +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC diagnostic pop +#endif + namespace { template @@ -235,6 +252,8 @@ namespace { using sample_types = ::testing::Types; TYPED_TEST_SUITE(fft_rt_test, sample_types); + // Already proved by the namespace-scope static_asserts above; this exists + // so the promise has a row in the test listing per profile. TYPED_TEST(fft_rt_test, TransformsAreNoexcept) { EXPECT_TRUE(transforms_are_noexcept()); } From 922491887d3b98211bb112bb93bc900fb1d4139d Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Thu, 17 Sep 2026 20:14:00 +0000 Subject: [PATCH 3/4] Register the parity targets through the bare-metal helper (#17 coupling) Rebased onto claude/wave1-stage1-legs. The two parity executables now go through tap_dsp_add_gtest_executable, so they are discovered hosted and one-shot under QEMU like tap_dsp_tests; the block no longer calls include(GoogleTest) or gtest_discover_tests itself. Consequences: - The 2^20 tests are compiled out (not GTEST_SKIP'ped) when TAP_DSP_PARITY_MAX_N is below 2^20, because the one-shot main counts a skip as a failed gate; the knob defaults to 4096 when CMAKE_CROSSCOMPILING and the toolchains' own cache value is honoured. - The oracle's closed-form sweep is capped by the same knob (tap_dsp_tests gets the definition), so no 65536-point buffers on the MPS2 legs. - The JUnit XML for the informational report is attached as a GTEST_OUTPUT ENVIRONMENT property via a ctest-time TEST_INCLUDE_FILES script, since the helper has no EXTRA_ARGS/PROPERTIES pass-through; hosted only. Verified on the host: GCC 13 and Clang 18 hosted (214/214, -Werror clean), -DTAP_DSP_BARE_METAL=ON host run (three one-shot markers, rc=0, selected=161/4/1, skipped=0, oracle and RT suites running under the negative filter), -DTAP_DSP_PARITY_MAX_N=4096 hosted (2^20 tests absent). Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- tests/CMakeLists.txt | 82 ++++++++++++++++++++------------- tests/test_fft_oracle.cpp | 10 +++- tests/test_fft_parity_ooura.cpp | 16 ++++--- 3 files changed, 67 insertions(+), 41 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 53e883e..93b23b3 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -145,12 +145,30 @@ target_link_libraries(tap_dsp_tests PRIVATE # BEGIN Stage 2a additions (docs/audit-fft-and-code-smells.md, Part 9): the # independent oracle and the real-time guard join the main executable; the # Ooura bit-identity gate gets its own targets so it can own its compiler flags. -# Keep this block self-contained and at the end of the file. +# Registration goes through tap_dsp_add_gtest_executable above, so the two +# targets are hosted (discovered) or bare-metal one-shot exactly like +# tap_dsp_tests. Keep this block self-contained and at the end of the file. # ============================================================================== +# Largest transform the FFT sweeps run. Hosted, the default covers every +# consumer geometry (2^20: five 8 MB buffers plus tables). Cross-compiled, the +# emulated Cortex-M legs (Part 10) get 4096, what fits the MPS2/MPS3 data +# regions; the toolchain files set the same value in the cache themselves and +# a plain cache set here never overrides theirs or a -D on the command line. +# The parity sweep, its 2^20 test (compiled out below the cap: a GTEST_SKIP is +# rc=1 on the target) and the oracle's closed-form sweep all read it. +if(CMAKE_CROSSCOMPILING) + set(_tap_dsp_parity_max_n_default 4096) +else() + set(_tap_dsp_parity_max_n_default 1048576) +endif() +set(TAP_DSP_PARITY_MAX_N ${_tap_dsp_parity_max_n_default} CACHE STRING + "Largest FFT size the Ooura parity and oracle sweeps run (power of two >= 4; 4096 on emulated targets)") + target_sources(tap_dsp_tests PRIVATE test_fft_oracle.cpp test_fft_rt.cpp) +target_compile_definitions(tap_dsp_tests PRIVATE TAP_DSP_PARITY_MAX_N=${TAP_DSP_PARITY_MAX_N}) # ------------------------------------------------------------------------------ # Two private reference builds of the vendored C. Neither is the tap_dsp_fft @@ -158,7 +176,8 @@ target_sources(tap_dsp_tests PRIVATE # TAP_DSP_FFT_CMSIS, where it also carries CMSIS objects), and the gate's whole # point is to compare against the C independently of how the library builds it. # The parity executables link one of these and NOT tap::dsp, so their -# rdft/rdft_f resolve here and no backend define (vDSP on Apple) reaches them. +# rdft/rdft_f resolve here and no backend define (vDSP on Apple, CMSIS on the +# M55) reaches them. # # Stage 2c note: fftsg.c moves to tests/reference/ooura/ (D6); this gate needs # fftsg_float.c to move with it, or its float side becomes port-vs-port. @@ -188,21 +207,11 @@ target_compile_options(tap_dsp_fft_reference_nocontract PRIVATE ${_tap_dsp_nocon add_library(tap_dsp_fft_reference_default STATIC ${_tap_dsp_reference_sources}) -# Largest transform the parity suite runs. Hosted, the default covers every -# consumer geometry (2^20: five 8 MB buffers plus tables). Cross-compiled, the -# emulated Cortex-M legs (Part 10) get 4096, the largest size that fits the -# MPS2/MPS3 data regions; a toolchain file may still set its own cache value, -# which this default does not override. -if(CMAKE_CROSSCOMPILING) - set(_tap_dsp_parity_max_n_default 4096) -else() - set(_tap_dsp_parity_max_n_default 1048576) -endif() -set(TAP_DSP_PARITY_MAX_N ${_tap_dsp_parity_max_n_default} CACHE STRING - "Largest FFT size the Ooura parity suite runs (power of two >= 4; 4096 on emulated targets)") - -# THE GATE: -ffp-contract=off on both sides, memcmp bit identity. -add_executable(tap_dsp_fft_parity test_fft_parity_ooura.cpp) +# THE GATE: -ffp-contract=off on both sides, memcmp bit identity. On the QEMU +# legs it runs whole (no MAIN_FILTER): at N <= 4096 the sweep is milliseconds. +tap_dsp_add_gtest_executable(tap_dsp_fft_parity + SOURCES test_fft_parity_ooura.cpp + LABELS parity) target_include_directories(tap_dsp_fft_parity PRIVATE ${PROJECT_SOURCE_DIR}/include) target_compile_features(tap_dsp_fft_parity PRIVATE cxx_std_20) target_compile_options(tap_dsp_fft_parity PRIVATE ${_tap_dsp_nocontract_cxx}) @@ -216,12 +225,13 @@ target_link_libraries(tap_dsp_fft_parity PRIVATE # assumed bound because a different fusion choice per butterfly stage # accumulates over log2 N stages; the number goes on the record per platform # instead (Part 6, item N3). Two channels carry it, because ctest hides the -# stdout of a passing test under --output-on-failure: (1) CI runs this label -# with -V as its own step (Part 13), (2) --gtest_output=xml writes the -# RecordProperty values to parity-ulp.xml in the build tree for upload. The -# report is ONE test covering both precisions so a single invocation produces -# the complete XML (each ctest entry is a separate process writing that file). -add_executable(tap_dsp_fft_parity_default_flags test_fft_parity_ooura.cpp) +# stdout of a passing test under --output-on-failure: (1) CI runs the parity +# label with -V as its own step (ci.yml, Part 13); (2) the JUnit XML below. +# The report is ONE test covering both precisions so a single invocation +# produces the complete XML. +tap_dsp_add_gtest_executable(tap_dsp_fft_parity_default_flags + SOURCES test_fft_parity_ooura.cpp + LABELS parity) target_include_directories(tap_dsp_fft_parity_default_flags PRIVATE ${PROJECT_SOURCE_DIR}/include) target_compile_features(tap_dsp_fft_parity_default_flags PRIVATE cxx_std_20) target_compile_definitions(tap_dsp_fft_parity_default_flags PRIVATE @@ -231,16 +241,22 @@ target_link_libraries(tap_dsp_fft_parity_default_flags PRIVATE tap_dsp_fft_reference_default tap_dsp_warnings) -include(GoogleTest) -target_link_libraries(tap_dsp_fft_parity PRIVATE GTest::gtest_main) -gtest_discover_tests(tap_dsp_fft_parity - DISCOVERY_TIMEOUT 120 - PROPERTIES TIMEOUT 900 LABELS parity) -target_link_libraries(tap_dsp_fft_parity_default_flags PRIVATE GTest::gtest_main) -gtest_discover_tests(tap_dsp_fft_parity_default_flags - EXTRA_ARGS --gtest_output=xml:${CMAKE_BINARY_DIR}/parity-ulp.xml - DISCOVERY_TIMEOUT 120 - PROPERTIES TIMEOUT 900 LABELS parity) +# JUnit XML for the informational report (hosted only: bare metal has no file +# system). gtest reads GTEST_OUTPUT from the environment at start-up, and the +# helper exposes no per-test properties or extra arguments, so the ENVIRONMENT +# property is attached the way gtest_discover_tests itself attaches its +# registrations: a script on the directory's TEST_INCLUDE_FILES, which ctest +# runs after the discovery script that defines the test (if(TEST ...) is not +# available in that context, so the call is unconditional: a discovery failure +# then errors here as well as in its own _NOT_BUILT entry). Should the helper +# grow an EXTRA_ARGS/PROPERTIES pass-through, this collapses to one line there. +if(NOT TAP_DSP_BARE_METAL) + set(_tap_dsp_parity_xml_script ${CMAKE_CURRENT_BINARY_DIR}/tap_dsp_fft_parity_default_flags_env.cmake) + file(WRITE ${_tap_dsp_parity_xml_script} + "set_tests_properties(fft_parity_ooura_default_flags.ReportsMaxUlpVersusOoura PROPERTIES\n" + " ENVIRONMENT \"GTEST_OUTPUT=xml:${CMAKE_BINARY_DIR}/parity-ulp.xml\")\n") + set_property(DIRECTORY APPEND PROPERTY TEST_INCLUDE_FILES ${_tap_dsp_parity_xml_script}) +endif() # ============================================================================== # END Stage 2a additions diff --git a/tests/test_fft_oracle.cpp b/tests/test_fft_oracle.cpp index e4b4fca..b37933c 100644 --- a/tests/test_fft_oracle.cpp +++ b/tests/test_fft_oracle.cpp @@ -53,6 +53,10 @@ #include "support/signals.h" #include "tap/dsp/fft.h" +#ifndef TAP_DSP_PARITY_MAX_N +#define TAP_DSP_PARITY_MAX_N (1 << 20) +#endif + namespace { // ------------------------------------------------------------------------ @@ -397,9 +401,13 @@ namespace { } } + /// The closed-form sweep: 4 .. min(65536, TAP_DSP_PARITY_MAX_N). The cap is + /// the same knob the parity gate uses (tests/CMakeLists.txt); the QEMU legs + /// set 4096 because a 65536-point double sweep needs several 512 KB + /// buffers that the MPS2 data region does not have. std::vector closed_form_sizes() { std::vector s; - for (std::size_t n = 4; n <= 65536; n *= 2) { + for (std::size_t n = 4; n <= 65536 && n <= static_cast(TAP_DSP_PARITY_MAX_N); n *= 2) { s.push_back(n); } return s; diff --git a/tests/test_fft_parity_ooura.cpp b/tests/test_fft_parity_ooura.cpp index 1dfd9eb..c6b6bf0 100644 --- a/tests/test_fft_parity_ooura.cpp +++ b/tests/test_fft_parity_ooura.cpp @@ -66,7 +66,8 @@ // // TAP_DSP_PARITY_MAX_N caps the sizes run, so the emulated QEMU legs can take // the suite at N <= 4096 (Part 9, Part 10) where 2^20 does not fit in RAM; -// the CMake block defaults it to 4096 when cross-compiling. +// the CMake block and the toolchain files default it to 4096 when +// cross-compiling, and the 2^20 test is compiled out (not skipped) below it. // // Stage 2c note. fftsg.c moves to tests/reference/ooura/ (Decision D6), and // the 2c text deletes fftsg_float.c. This gate needs BOTH files: the float @@ -290,7 +291,8 @@ namespace { return sizes; } - constexpr std::size_t k_large_n = std::size_t{1} << 20; + // Used by the gate only when it is compiled in (see the #if below). + [[maybe_unused]] constexpr std::size_t k_large_n = std::size_t{1} << 20; #if !defined(TAP_DSP_PARITY_INFORMATIONAL) @@ -338,13 +340,12 @@ namespace { // One run at 2^20 — the largest geometry any consumer uses (AmbiTap's // long-partition convolution) — kept as its own test so its runtime shows - // separately and it can be excluded on emulated targets by name or by - // TAP_DSP_PARITY_MAX_N. + // separately. COMPILED OUT, not skipped, when TAP_DSP_PARITY_MAX_N is + // below 2^20: the bare-metal one-shot main counts a GTEST_SKIP as a failed + // gate (tests/bare_metal_main.cpp), and the QEMU legs run at 4096. +#if TAP_DSP_PARITY_MAX_N >= (1 << 20) template void expect_identical_large() { - if (k_large_n > static_cast(TAP_DSP_PARITY_MAX_N)) { - GTEST_SKIP() << "2^20 exceeds TAP_DSP_PARITY_MAX_N=" << TAP_DSP_PARITY_MAX_N; - } for (const material m : k_materials) { expect_identical(k_large_n, m, true); expect_identical(k_large_n, m, false); @@ -361,6 +362,7 @@ namespace { TEST(fft_parity_ooura, LargeTransformIsBitIdenticalToOouraFloat) { expect_identical_large(); } +#endif // TAP_DSP_PARITY_MAX_N >= (1 << 20) #else // TAP_DSP_PARITY_INFORMATIONAL From b303e706e856722651e9290b478fcb7473c9cff4 Mon Sep 17 00:00:00 2001 From: Timothy Place Date: Thu, 17 Sep 2026 20:25:30 +0000 Subject: [PATCH 4/4] Confine the float oracle sweeps to the CMSIS backend's 32..4096 on the M55 The first QEMU run of the rebased suite HardFaulted in fft_oracle_test/0.ImpulseIsFlatAtEverySize at N = 4 on the cortex-m55 leg, the one leg where basic_real_fft routes through CMSIS-DSP: arm_rfft_fast_init_f32 accepts only 32..4096 and fft.h's wrapper does not check its return status, so fft.h's ">= 4" contract is undefined behaviour there below 32 (the existing battery starts at 64, which is why it never showed). The oracle's profile now carries the engine's size range (k_min_n / k_max_n; 32..4096 for float under TAP_DSP_FFT_CMSIS, 4..2^20 otherwise) and every sweep is built from it. The wrapper defect itself is recorded for the fft.h owner, not papered over here. Also: the informational table printed "zu" for N on newlib-nano; %lu now. Hosted GCC 13 / Clang 18: 214/214, -Werror clean; host one-shot mode: three markers rc=0, selected=161/4/1, skipped=0. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_019ZPTzNxo5Fe4EtpXXKf7Sy --- tests/test_fft_oracle.cpp | 115 ++++++++++++++++++++------------ tests/test_fft_parity_ooura.cpp | 6 +- 2 files changed, 78 insertions(+), 43 deletions(-) diff --git a/tests/test_fft_oracle.cpp b/tests/test_fft_oracle.cpp index b37933c..76845d8 100644 --- a/tests/test_fft_oracle.cpp +++ b/tests/test_fft_oracle.cpp @@ -125,26 +125,51 @@ namespace { template struct profile; + // Size range a profile's engine accepts. fft.h promises any power of two + // >= 4 for Ooura, and the double profile is always Ooura. The float + // profile is whatever backend the build selected: under TAP_DSP_FFT_CMSIS + // (the M55 leg) CMSIS-DSP's arm_rfft_fast_init_f32 accepts 32..4096 only + // (arm_rfft_fast_init_f32.c, the switch at the end) and fft.h's wrapper + // does not check its return status, so a size outside that range is + // undefined behaviour — a HardFault at N = 4 on the QEMU M55 leg is how + // this was found. Until fft.h rejects or falls back on those sizes (a + // finding for the fft.h owner, not this file), the float sweeps on that + // backend run over the range the backend supports. vDSP (macOS) takes the + // full range. + constexpr std::size_t k_ooura_min_n = 4; + constexpr std::size_t k_ooura_max_n = std::size_t{1} << 20; +#if defined(TAP_DSP_FFT_CMSIS) + constexpr std::size_t k_float_backend_min_n = 32; + constexpr std::size_t k_float_backend_max_n = 4096; +#else + constexpr std::size_t k_float_backend_min_n = k_ooura_min_n; + constexpr std::size_t k_float_backend_max_n = k_ooura_max_n; +#endif + template <> struct profile { - static constexpr double k_epsilon = std::numeric_limits::epsilon(); - static constexpr double k_full_scale = 1.0; ///< closed-form drive amplitude - static double to_double(double v) { return v; } - static double from_double(double v) { return v; } - static double forward_scale(std::size_t) { return 1.0; } ///< engine forward = scale * DFT - static double inverse_scale(std::size_t) { return 1.0; } ///< engine inverse = scale * unnormalized - static double tolerance(std::size_t n, double norm2) { return higham_tolerance(k_epsilon, n, norm2); } + static constexpr double k_epsilon = std::numeric_limits::epsilon(); + static constexpr double k_full_scale = 1.0; ///< closed-form drive amplitude + static constexpr std::size_t k_min_n = k_ooura_min_n; + static constexpr std::size_t k_max_n = k_ooura_max_n; + static double to_double(double v) { return v; } + static double from_double(double v) { return v; } + static double forward_scale(std::size_t) { return 1.0; } ///< engine forward = scale * DFT + static double inverse_scale(std::size_t) { return 1.0; } ///< engine inverse = scale * unnormalized + static double tolerance(std::size_t n, double norm2) { return higham_tolerance(k_epsilon, n, norm2); } }; template <> struct profile { - static constexpr double k_epsilon = std::numeric_limits::epsilon(); - static constexpr double k_full_scale = 1.0; - static double to_double(float v) { return static_cast(v); } - static float from_double(double v) { return static_cast(v); } - static double forward_scale(std::size_t) { return 1.0; } - static double inverse_scale(std::size_t) { return 1.0; } - static double tolerance(std::size_t n, double norm2) { return higham_tolerance(k_epsilon, n, norm2); } + static constexpr double k_epsilon = std::numeric_limits::epsilon(); + static constexpr double k_full_scale = 1.0; + static constexpr std::size_t k_min_n = k_float_backend_min_n; + static constexpr std::size_t k_max_n = k_float_backend_max_n; + static double to_double(float v) { return static_cast(v); } + static float from_double(double v) { return static_cast(v); } + static double forward_scale(std::size_t) { return 1.0; } + static double inverse_scale(std::size_t) { return 1.0; } + static double tolerance(std::size_t n, double norm2) { return higham_tolerance(k_epsilon, n, norm2); } }; using oracle_types = ::testing::Types; @@ -401,24 +426,32 @@ namespace { } } - /// The closed-form sweep: 4 .. min(65536, TAP_DSP_PARITY_MAX_N). The cap is - /// the same knob the parity gate uses (tests/CMakeLists.txt); the QEMU legs - /// set 4096 because a 65536-point double sweep needs several 512 KB - /// buffers that the MPS2 data region does not have. - std::vector closed_form_sizes() { + /// Powers of two from the profile's minimum up to min(limit, + /// TAP_DSP_PARITY_MAX_N, the profile's maximum). The cap is the same knob + /// the parity gate uses (tests/CMakeLists.txt); the QEMU legs set 4096 + /// because a 65536-point double sweep needs several 512 KB buffers that + /// the MPS2 data region does not have. + template + std::vector sizes_up_to(std::size_t limit) { + const std::size_t top = + std::min({limit, static_cast(TAP_DSP_PARITY_MAX_N), profile::k_max_n}); std::vector s; - for (std::size_t n = 4; n <= 65536 && n <= static_cast(TAP_DSP_PARITY_MAX_N); n *= 2) { + for (std::size_t n = profile::k_min_n; n <= top; n *= 2) { s.push_back(n); } return s; } + /// The closed-form sweep: up to 65536. + template + std::vector closed_form_sizes() { + return sizes_up_to(65536); + } + + /// The compensated-DFT sizes: up to 256 (see "What the bound cannot see"). + template std::vector dft_sizes() { - std::vector s; - for (std::size_t n = 4; n <= 256; n *= 2) { - s.push_back(n); - } - return s; + return sizes_up_to(256); } /// Bin used for the on-bin materials: n/8, or 1 below n = 8. Always in @@ -472,7 +505,7 @@ namespace { TYPED_TEST(fft_oracle_test, ImpulseIsFlatAtEverySize) { const double a = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { std::vector x(n, 0.0); x[0] = a; std::vector expected(n, 0.0); @@ -487,7 +520,7 @@ namespace { TYPED_TEST(fft_oracle_test, DcLandsAsNInSlotZero) { const double a = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { std::vector x(n, a); std::vector expected(n, 0.0); expected[0] = a * static_cast(n); @@ -497,7 +530,7 @@ namespace { TYPED_TEST(fft_oracle_test, NyquistAlternationLandsAsNInSlotOne) { const double a = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { std::vector x(n); for (std::size_t j = 0; j < n; ++j) { x[j] = (j % 2 == 0) ? a : -a; @@ -510,7 +543,7 @@ namespace { TYPED_TEST(fft_oracle_test, OnBinCosineIsPlusHalfNInTheRealSlot) { const double a = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { const std::size_t k = tone_bin(n); const auto x = tap::dsp::test::tone(n, static_cast(k), a, 0.0); std::vector expected(n, 0.0); @@ -523,7 +556,7 @@ namespace { // (it would be -N/2 in the engineering convention exp(-2*pi*i/N)). TYPED_TEST(fft_oracle_test, OnBinSineIsPlusHalfNInTheImaginarySlot) { const double a = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { const std::size_t k = tone_bin(n); // sin(w j) = cos(w j - pi/2) const auto x = tap::dsp::test::tone(n, static_cast(k), a, -std::numbers::pi / 2.0); @@ -535,7 +568,7 @@ namespace { TYPED_TEST(fft_oracle_test, TwoToneSuperposesLinearly) { const double a = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { if (n < 8) { continue; // needs two distinct complex bins } @@ -562,7 +595,7 @@ namespace { TYPED_TEST(fft_oracle_test, InverseOfFlatSpectrumIsHalfNImpulse) { const double amp = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { std::vector a(n, 0.0); a[0] = amp; a[1] = amp; @@ -577,7 +610,7 @@ namespace { TYPED_TEST(fft_oracle_test, InverseOfDcSpectrumIsHalfNConstant) { const double amp = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { std::vector a(n, 0.0); a[0] = amp; std::vector expected(n, 0.5 * amp); @@ -587,7 +620,7 @@ namespace { TYPED_TEST(fft_oracle_test, InverseOfNyquistSpectrumIsHalfNAlternation) { const double amp = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { std::vector a(n, 0.0); a[1] = amp; std::vector expected(n); @@ -600,7 +633,7 @@ namespace { TYPED_TEST(fft_oracle_test, InverseOfOnBinSpectrumIsHalfNTone) { const double amp = profile::k_full_scale; - for (const std::size_t n : closed_form_sizes()) { + for (const std::size_t n : closed_form_sizes()) { const std::size_t k = tone_bin(n); std::vector a(n, 0.0); a[2 * k] = amp; // cosine part @@ -636,7 +669,7 @@ namespace { // multiply (exact, power of two), so ~1.5 eps relative on unit-scale data; // constant 8 eps; measured maximum 1.1e-16 = 0.5 eps. TEST(fft_oracle_self_check, TwiddlesMatchLibmToDoubleRounding) { - for (const std::size_t n : dft_sizes()) { + for (const std::size_t n : dft_sizes()) { const compensated_dft oracle(n); for (std::size_t m = 0; m < n; ++m) { const double theta = 2.0 * std::numbers::pi * static_cast(m) / static_cast(n); @@ -653,7 +686,7 @@ namespace { } TEST(fft_oracle_self_check, RoundTripIsIdentityToDoubleRounding) { - for (const std::size_t n : dft_sizes()) { + for (const std::size_t n : dft_sizes()) { const compensated_dft oracle(n); const std::vector x = tap::dsp::test::random_signal(n, 0xC0FFEEu); const std::vector a = oracle.forward(x); @@ -681,7 +714,7 @@ namespace { } TYPED_TEST(fft_oracle_test, ForwardMatchesCompensatedDftOnBroadband) { - for (const std::size_t n : dft_sizes()) { + for (const std::size_t n : dft_sizes()) { // Drawn in the profile first so the oracle sees exactly the bits // the engine sees; from_doubles is then the identity. const auto x = to_doubles(tap::dsp::test::random_signal(n, 0x2545F491u)); @@ -690,7 +723,7 @@ namespace { } TYPED_TEST(fft_oracle_test, ForwardMatchesCompensatedDftOnOffBinTone) { - for (const std::size_t n : dft_sizes()) { + for (const std::size_t n : dft_sizes()) { // Off-bin (k + 0.37) so every bin carries leakage and none is // numerically empty: the opposite regime from the closed forms. const auto x = @@ -700,7 +733,7 @@ namespace { } TYPED_TEST(fft_oracle_test, InverseMatchesCompensatedDftOnBroadbandSpectrum) { - for (const std::size_t n : dft_sizes()) { + for (const std::size_t n : dft_sizes()) { // An arbitrary packed spectrum, not one produced by a forward, so // the inverse is judged on its own definition. const auto a = to_doubles(tap::dsp::test::random_signal(n, 0x1D872B41u)); @@ -709,7 +742,7 @@ namespace { } TYPED_TEST(fft_oracle_test, InverseMatchesCompensatedDftOnToneSpectrum) { - for (const std::size_t n : dft_sizes()) { + for (const std::size_t n : dft_sizes()) { const compensated_dft oracle(n); const auto x = to_doubles(tap::dsp::test::tone(n, static_cast(tone_bin(n)) + 0.37, 0.8, 1.1)); diff --git a/tests/test_fft_parity_ooura.cpp b/tests/test_fft_parity_ooura.cpp index c6b6bf0..410f11e 100644 --- a/tests/test_fft_parity_ooura.cpp +++ b/tests/test_fft_parity_ooura.cpp @@ -401,8 +401,10 @@ namespace { } } overall = std::max({overall, fwd, inv}); - std::printf("[ ulp ] %10zu %14llu %14llu %s / %s\n", n, static_cast(fwd), - static_cast(inv), fwd_worst, inv_worst); + // %lu, not %zu: newlib-nano's printf on the QEMU legs prints "zu". + std::printf("[ ulp ] %10lu %14llu %14llu %s / %s\n", static_cast(n), + static_cast(fwd), static_cast(inv), fwd_worst, + inv_worst); ::testing::Test::RecordProperty(std::string(precision_name()) + "_N" + std::to_string(n) + "_fwd", std::to_string(fwd)); ::testing::Test::RecordProperty(std::string(precision_name()) + "_N" + std::to_string(n) + "_inv",