Merge pull request #43 from falseywinchnet/c93rxr-main/implement-avx2-pathway-for-odd-fft-transform

Add AVX2 forward-combine path to BODFT kernel

Author: falseywinchnet

Makefile

@@ -66,10 +66,11 @@ SFDR_PROBE := $(BUILD_DIR)/tests/bfft_fftw_sfdr_probe
 BH7_PROBE := $(BUILD_DIR)/tests/bfft_fftw_sfdr_bh7_probe
 LIBRARY_COMPARE_PROBE := $(BUILD_DIR)/tests/bfft_library_compare_probe
 AVX2_ASM := $(BUILD_DIR)/src/bfft_avx2.s
+BODFT_AVX2_ASM := $(BUILD_DIR)/src/bodft_avx2.s
 SSE2_ASM := $(BUILD_DIR)/src/bfft_sse2.s
 ASM_OUTPUTS :=
 ifneq ($(AVX2_FLAGS),)
-  ASM_OUTPUTS += $(AVX2_ASM)
+  ASM_OUTPUTS += $(AVX2_ASM) $(BODFT_AVX2_ASM)
 endif
 ifneq ($(SSE2_FLAGS),)
   ASM_OUTPUTS += $(SSE2_ASM)
@@ -109,6 +110,9 @@ asm-check: $(ASM_OUTPUTS)
 $(AVX2_ASM): $(SRC) include/bfft/bfft.h src/detail/bruun_kernel.hpp | $(BUILD_DIR)
 	$(CXX) $(LIB_CPPFLAGS) $(CXXFLAGS) $(AVX2_FLAGS) -S -fverbose-asm $< -o $@
 
+$(BODFT_AVX2_ASM): $(BODFT_SRC) include/bfft/bodft.h include/bfft/bfft.h src/detail/bodft_kernel.hpp src/detail/bruun_kernel.hpp | $(BUILD_DIR)
+	$(CXX) $(LIB_CPPFLAGS) $(CXXFLAGS) $(AVX2_FLAGS) -S -fverbose-asm $< -o $@
+
 $(SSE2_ASM): $(SRC) include/bfft/bfft.h src/detail/bruun_kernel.hpp | $(BUILD_DIR)
 	$(CXX) $(LIB_CPPFLAGS) $(CXXFLAGS) $(SSE2_FLAGS) -S -fverbose-asm $< -o $@
 

src/detail/bodft_avx2_notes.md

@@ -0,0 +1,23 @@
+# BODFT AVX2 pathway notes
+
+## 2026-06-23
+
+- Added a 256-bit AVX2/FMA double-precision forward combine path for the odd-frequency transform.
+- The AVX2 code handles four paired radix-4 positions per iteration, twice the double-lane width of the SSE2/NEON path.
+- The path uses structure-of-arrays lane registers for child spectra and twiddle tables, performs the three child twiddle complex multiplies with FMA instructions, then stores the forward and conjugate-partner outputs through one packing helper.
+- Kept the existing 128-bit SIMD double path as the fallback tail and non-AVX2 implementation. Float remains on the portable 128-bit SIMD path because its current path already handles four positions per iteration and needs a separate eight-position AVX2 packing pass to be worthwhile.
+- Added a Makefile assembly target for `src/bodft.cpp` so `make asm-check` emits both the regular BFFT AVX2 assembly and the BODFT AVX2 assembly.
+
+## Validation notes
+
+- Correctness was checked with `make test BUILD_DIR=build-avx2 CXXWARNFLAGS='-Wall -Wextra -Wpedantic -Werror'`.
+- BODFT-specific correctness and timing were checked with `make build-avx2/examples/bodft_benchmark BUILD_DIR=build-avx2 && build-avx2/examples/bodft_benchmark`.
+- Assembly was checked with `make asm-check BUILD_DIR=build-avx2` and by scanning `build-avx2/src/bodft_avx2.s` for AVX/FMA instructions.
+
+## Performance observations from this container
+
+- Backend reported by the benchmark: `avx2-fma-256`.
+- Double forward timings from the AVX2 path were approximately 0.83 us at N=256, 3.88 us at N=1024, 19.2 us at N=4096, 458 us at N=65536, and 18.3 ms at N=1048576.
+- Float forward timings remained on the existing 128-bit path and were approximately 0.71 us at N=256, 3.13 us at N=1024, 14.7 us at N=4096, 313 us at N=65536, and 11.9 ms at N=1048576.
+- The next performance task is an eight-position AVX2 float combine and, after that, an AVX2 inverse combine if inverse parity becomes a requirement.
+- For a same-container point comparison, `build-avx2/examples/benchmark 65536 200` measured the regular BFFT native double path at about 217 us for N=65536; the BODFT double forward path measured about 458 us for N=65536 in the broader BODFT benchmark. The gap is now concentrated in BODFT's transform-specific data movement/scratch schedule and unvectorized inverse rather than missing AVX2 forward-combine code.

src/detail/bodft_kernel.hpp

@@ -90,6 +90,81 @@ static inline CT cmuli(CT a) { return CT{-a.im, a.re}; }
 #  endif
 #endif
 
+
+#if BRUUN_LEVEL >= 2
+static inline void bodft_uzp4_pd(const complex_t* RESTRICT src, int k,
+                                 __m256d& re, __m256d& im) {
+    const __m256d a = _mm256_loadu_pd(&src[k].re);
+    const __m256d b = _mm256_loadu_pd(&src[k + 2].re);
+    const __m256d lo = _mm256_unpacklo_pd(a, b);
+    const __m256d hi = _mm256_unpackhi_pd(a, b);
+    re = _mm256_permute4x64_pd(lo, 0xD8);
+    im = _mm256_permute4x64_pd(hi, 0xD8);
+}
+
+static inline void bodft_store4_pd(complex_t* RESTRICT dst, int k,
+                                   __m256d re, __m256d im) {
+    const __m256d lo = _mm256_unpacklo_pd(re, im);
+    const __m256d hi = _mm256_unpackhi_pd(re, im);
+    _mm256_storeu_pd(&dst[k].re,     _mm256_permute2f128_pd(lo, hi, 0x20));
+    _mm256_storeu_pd(&dst[k + 2].re, _mm256_permute2f128_pd(lo, hi, 0x31));
+}
+
+static inline __m256d bodft_rev4_pd(__m256d v) {
+    return _mm256_permute4x64_pd(v, 0x1B);
+}
+
+static inline int bodft_combine_fwd_avx2_f64(const complex_t* RESTRICT tab,
+                                             const complex_t* RESTRICT tab2,
+                                             const complex_t* RESTRICT tab3,
+                                             const complex_t* RESTRICT c0,
+                                             const complex_t* RESTRICT c1,
+                                             const complex_t* RESTRICT c2,
+                                             const complex_t* RESTRICT c3,
+                                             complex_t* RESTRICT out,
+                                             int M, int half) {
+    const __m256d zero = _mm256_setzero_pd();
+    int k = 0;
+    for (; k + 4 <= half; k += 4) {
+        __m256d tre, tim, t2re, t2im, t3re, t3im;
+        __m256d c0re, c0im, c1re, c1im, c2re, c2im, c3re, c3im;
+        bodft_uzp4_pd(tab,  k, tre, tim);
+        bodft_uzp4_pd(tab2, k, t2re, t2im);
+        bodft_uzp4_pd(tab3, k, t3re, t3im);
+        bodft_uzp4_pd(c0, k, c0re, c0im);
+        bodft_uzp4_pd(c1, k, c1re, c1im);
+        bodft_uzp4_pd(c2, k, c2re, c2im);
+        bodft_uzp4_pd(c3, k, c3re, c3im);
+
+        const __m256d b0re = c0re, b0im = c0im;
+        const __m256d b1re = _mm256_fmsub_pd(c1re, tre, _mm256_mul_pd(c1im, tim));
+        const __m256d b1im = _mm256_fmadd_pd(c1re, tim, _mm256_mul_pd(c1im, tre));
+        const __m256d b2re = _mm256_fmsub_pd(c2re, t2re, _mm256_mul_pd(c2im, t2im));
+        const __m256d b2im = _mm256_fmadd_pd(c2re, t2im, _mm256_mul_pd(c2im, t2re));
+        const __m256d b3re = _mm256_fmsub_pd(c3re, t3re, _mm256_mul_pd(c3im, t3im));
+        const __m256d b3im = _mm256_fmadd_pd(c3re, t3im, _mm256_mul_pd(c3im, t3re));
+
+        const __m256d e0re = _mm256_add_pd(b0re, b2re), e0im = _mm256_add_pd(b0im, b2im);
+        const __m256d e1re = _mm256_sub_pd(b0re, b2re), e1im = _mm256_sub_pd(b0im, b2im);
+        const __m256d o0re = _mm256_add_pd(b1re, b3re), o0im = _mm256_add_pd(b1im, b3im);
+        const __m256d o1re = _mm256_sub_pd(b1re, b3re), o1im = _mm256_sub_pd(b1im, b3im);
+
+        const __m256d yk_re = _mm256_add_pd(e0re, o0re), yk_im = _mm256_add_pd(e0im, o0im);
+        const __m256d ykM_re = _mm256_add_pd(e1re, o1im), ykM_im = _mm256_sub_pd(e1im, o1re);
+        const __m256d ykp_re = _mm256_sub_pd(e1re, o1im);
+        const __m256d ykp_im = _mm256_sub_pd(zero, _mm256_add_pd(e1im, o1re));
+        const __m256d ykpM_re = _mm256_sub_pd(e0re, o0re);
+        const __m256d ykpM_im = _mm256_sub_pd(o0im, e0im);
+
+        bodft_store4_pd(out, k, yk_re, yk_im);
+        bodft_store4_pd(out, k + M, ykM_re, ykM_im);
+        bodft_store4_pd(out, M - 4 - k, bodft_rev4_pd(ykp_re), bodft_rev4_pd(ykp_im));
+        bodft_store4_pd(out, M - 4 - k + M, bodft_rev4_pd(ykpM_re), bodft_rev4_pd(ykpM_im));
+    }
+    return k;
+}
+#endif
+
 // ---------------------------------------------------------------------------
 // Forward double combine, 128-bit SoA SIMD specialization. Two output positions
 // k, k+1 per iteration; children and twiddles transpose into real/imag lane
@@ -234,9 +309,22 @@ static inline void combine_fwd(const CT* RESTRICT tab, const CT* RESTRICT tab2,
                                const CT* RESTRICT c3, CT* RESTRICT out,
                                int M, int half) {
     int k = 0;
+#if BRUUN_LEVEL >= 2
+    if constexpr (sizeof(RT) == 8) {
+        k = bodft_combine_fwd_avx2_f64(
+            reinterpret_cast<const complex_t*>(tab),
+            reinterpret_cast<const complex_t*>(tab2),
+            reinterpret_cast<const complex_t*>(tab3),
+            reinterpret_cast<const complex_t*>(c0),
+            reinterpret_cast<const complex_t*>(c1),
+            reinterpret_cast<const complex_t*>(c2),
+            reinterpret_cast<const complex_t*>(c3),
+            reinterpret_cast<complex_t*>(out), M, half);
+    }
+#endif
 #if BRUUN_LEVEL >= 1
     if constexpr (sizeof(RT) == 8) {
-        k = bodft_combine_fwd_simd_f64(
+        if (k == 0) k = bodft_combine_fwd_simd_f64(
             reinterpret_cast<const complex_t*>(tab),
             reinterpret_cast<const complex_t*>(tab2),
             reinterpret_cast<const complex_t*>(tab3),