Merge pull request #72 from falseywinchnet/yxlmq0-main/optimize-performance-with-avx2/sse2

Add AVX2 float BODFT forward combine (8-position, FMA)

Author: falseywinchnet

src/detail/bodft_avx2_notes.md

@@ -29,3 +29,11 @@
 - Tightened the forward scratch schedule so the final forward combine writes directly into the caller output buffer instead of ping-ponging once more and copying `N/2` complex values at the end.
 - BODFT double inverse timing improved in this container from about 748 us to about 611 us at N=65536, and from about 22.7 ms to about 20.6 ms at N=1048576. Smaller sizes saw larger relative gains where the final-copy removal and inverse AVX2 combine reduce overhead together.
 - The remaining gap to the regular BFFT native path is now mostly scratch/data-order locality and float-specific AVX2 work rather than the absence of a double inverse vector combine.
+
+## Float forward combine pass
+
+- Added an eight-position AVX2/FMA single-precision forward combine path for the odd-frequency transform.
+- The float path now deinterleaves eight packed complex twiddles/children into 256-bit SoA real/imaginary lanes, does all three child twiddle products with `vfmaddps`/`vfmsubps`, and writes the four paired output blocks through an AVX2 interleave helper.
+- Partner bins are lane-reversed in registers with an AVX2 permute before storing, preserving the scalar/128-bit conjugate-pair layout while doubling the 128-bit float combine width on AVX2 hosts.
+- The existing 128-bit SSE2/NEON float combine remains the fallback for non-AVX2 builds and for tiny levels that cannot fill an eight-position vector.
+- Correctness was checked with the regular CTest suite and the standalone BODFT benchmark on an AVX2/FMA build; a separate SSE2/no-AVX build was also checked to keep the fallback path intact.

src/detail/bodft_kernel.hpp

@@ -164,6 +164,81 @@ static inline int bodft_combine_fwd_avx2_f64(const complex_t* RESTRICT tab,
     return k;
 }
 
+
+static inline void bodft_uzp8_ps(const complex_f32_t* RESTRICT src, int k,
+                                 __m256& re, __m256& im) {
+    const __m256 a = _mm256_loadu_ps(&src[k].re);
+    const __m256 b = _mm256_loadu_ps(&src[k + 4].re);
+    const __m256 lo = _mm256_unpacklo_ps(a, b);
+    const __m256 hi = _mm256_unpackhi_ps(a, b);
+    const __m256 re_pairs = _mm256_shuffle_ps(lo, hi, _MM_SHUFFLE(1, 0, 1, 0));
+    const __m256 im_pairs = _mm256_shuffle_ps(lo, hi, _MM_SHUFFLE(3, 2, 3, 2));
+    re = _mm256_permutevar8x32_ps(re_pairs, _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7));
+    im = _mm256_permutevar8x32_ps(im_pairs, _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7));
+}
+
+static inline void bodft_store8_ps(complex_f32_t* RESTRICT dst, int k,
+                                   __m256 re, __m256 im) {
+    const __m256 lo = _mm256_unpacklo_ps(re, im);
+    const __m256 hi = _mm256_unpackhi_ps(re, im);
+    _mm256_storeu_ps(&dst[k].re,     _mm256_permute2f128_ps(lo, hi, 0x20));
+    _mm256_storeu_ps(&dst[k + 4].re, _mm256_permute2f128_ps(lo, hi, 0x31));
+}
+
+static inline __m256 bodft_rev8_ps(__m256 v) {
+    return _mm256_permutevar8x32_ps(v, _mm256_setr_epi32(7, 6, 5, 4, 3, 2, 1, 0));
+}
+
+static inline int bodft_combine_fwd_avx2_f32(const complex_f32_t* RESTRICT tab,
+                                             const complex_f32_t* RESTRICT tab2,
+                                             const complex_f32_t* RESTRICT tab3,
+                                             const complex_f32_t* RESTRICT c0,
+                                             const complex_f32_t* RESTRICT c1,
+                                             const complex_f32_t* RESTRICT c2,
+                                             const complex_f32_t* RESTRICT c3,
+                                             complex_f32_t* RESTRICT out,
+                                             int M, int half) {
+    const __m256 zero = _mm256_setzero_ps();
+    int k = 0;
+    for (; k + 8 <= half; k += 8) {
+        __m256 tre, tim, t2re, t2im, t3re, t3im;
+        __m256 c0re, c0im, c1re, c1im, c2re, c2im, c3re, c3im;
+        bodft_uzp8_ps(tab,  k, tre, tim);
+        bodft_uzp8_ps(tab2, k, t2re, t2im);
+        bodft_uzp8_ps(tab3, k, t3re, t3im);
+        bodft_uzp8_ps(c0, k, c0re, c0im);
+        bodft_uzp8_ps(c1, k, c1re, c1im);
+        bodft_uzp8_ps(c2, k, c2re, c2im);
+        bodft_uzp8_ps(c3, k, c3re, c3im);
+
+        const __m256 b0re = c0re, b0im = c0im;
+        const __m256 b1re = _mm256_fmsub_ps(c1re, tre, _mm256_mul_ps(c1im, tim));
+        const __m256 b1im = _mm256_fmadd_ps(c1re, tim, _mm256_mul_ps(c1im, tre));
+        const __m256 b2re = _mm256_fmsub_ps(c2re, t2re, _mm256_mul_ps(c2im, t2im));
+        const __m256 b2im = _mm256_fmadd_ps(c2re, t2im, _mm256_mul_ps(c2im, t2re));
+        const __m256 b3re = _mm256_fmsub_ps(c3re, t3re, _mm256_mul_ps(c3im, t3im));
+        const __m256 b3im = _mm256_fmadd_ps(c3re, t3im, _mm256_mul_ps(c3im, t3re));
+
+        const __m256 e0re = _mm256_add_ps(b0re, b2re), e0im = _mm256_add_ps(b0im, b2im);
+        const __m256 e1re = _mm256_sub_ps(b0re, b2re), e1im = _mm256_sub_ps(b0im, b2im);
+        const __m256 o0re = _mm256_add_ps(b1re, b3re), o0im = _mm256_add_ps(b1im, b3im);
+        const __m256 o1re = _mm256_sub_ps(b1re, b3re), o1im = _mm256_sub_ps(b1im, b3im);
+
+        const __m256 yk_re = _mm256_add_ps(e0re, o0re), yk_im = _mm256_add_ps(e0im, o0im);
+        const __m256 ykM_re = _mm256_add_ps(e1re, o1im), ykM_im = _mm256_sub_ps(e1im, o1re);
+        const __m256 ykp_re = _mm256_sub_ps(e1re, o1im);
+        const __m256 ykp_im = _mm256_sub_ps(zero, _mm256_add_ps(e1im, o1re));
+        const __m256 ykpM_re = _mm256_sub_ps(e0re, o0re);
+        const __m256 ykpM_im = _mm256_sub_ps(o0im, e0im);
+
+        bodft_store8_ps(out, k, yk_re, yk_im);
+        bodft_store8_ps(out, k + M, ykM_re, ykM_im);
+        bodft_store8_ps(out, M - 8 - k, bodft_rev8_ps(ykp_re), bodft_rev8_ps(ykp_im));
+        bodft_store8_ps(out, M - 8 - k + M, bodft_rev8_ps(ykpM_re), bodft_rev8_ps(ykpM_im));
+    }
+    return k;
+}
+
 static inline int bodft_combine_inv_avx2_f64(const complex_t* RESTRICT tab,
                                              const complex_t* RESTRICT tab2,
                                              const complex_t* RESTRICT tab3,
@@ -384,6 +459,19 @@ static inline void combine_fwd(const CT* RESTRICT tab, const CT* RESTRICT tab2,
             reinterpret_cast<complex_t*>(out), M, half);
     }
 #endif
+#if BRUUN_LEVEL >= 2
+    if constexpr (sizeof(RT) == 4) {
+        k = bodft_combine_fwd_avx2_f32(
+            reinterpret_cast<const complex_f32_t*>(tab),
+            reinterpret_cast<const complex_f32_t*>(tab2),
+            reinterpret_cast<const complex_f32_t*>(tab3),
+            reinterpret_cast<const complex_f32_t*>(c0),
+            reinterpret_cast<const complex_f32_t*>(c1),
+            reinterpret_cast<const complex_f32_t*>(c2),
+            reinterpret_cast<const complex_f32_t*>(c3),
+            reinterpret_cast<complex_f32_t*>(out), M, half);
+    }
+#endif
 #if BRUUN_LEVEL >= 1
     if constexpr (sizeof(RT) == 8) {
         if (k == 0) k = bodft_combine_fwd_simd_f64(
@@ -396,15 +484,17 @@ static inline void combine_fwd(const CT* RESTRICT tab, const CT* RESTRICT tab2,
             reinterpret_cast<const complex_t*>(c3),
             reinterpret_cast<complex_t*>(out), M, half);
     } else {
-        k = bodft_combine_fwd_simd_f32(
-            reinterpret_cast<const complex_f32_t*>(tab),
-            reinterpret_cast<const complex_f32_t*>(tab2),
-            reinterpret_cast<const complex_f32_t*>(tab3),
-            reinterpret_cast<const complex_f32_t*>(c0),
-            reinterpret_cast<const complex_f32_t*>(c1),
-            reinterpret_cast<const complex_f32_t*>(c2),
-            reinterpret_cast<const complex_f32_t*>(c3),
-            reinterpret_cast<complex_f32_t*>(out), M, half);
+        if (k == 0) {
+            k = bodft_combine_fwd_simd_f32(
+                reinterpret_cast<const complex_f32_t*>(tab),
+                reinterpret_cast<const complex_f32_t*>(tab2),
+                reinterpret_cast<const complex_f32_t*>(tab3),
+                reinterpret_cast<const complex_f32_t*>(c0),
+                reinterpret_cast<const complex_f32_t*>(c1),
+                reinterpret_cast<const complex_f32_t*>(c2),
+                reinterpret_cast<const complex_f32_t*>(c3),
+                reinterpret_cast<complex_f32_t*>(out), M, half);
+        }
     }
 #endif
     for (; k < half; ++k) {