Skip to content

Commit cfb88b3

Browse files
Merge pull request #73 from falseywinchnet/fwysqa-main/optimize-performance-with-avx2/sse2
AVX2: eight-wide f32 BODFT forward/inverse combines, FMA intrinsics and SIMD tail scheduling improvements
2 parents aa658f9 + 356ea25 commit cfb88b3

5 files changed

Lines changed: 358 additions & 45 deletions

File tree

src/detail/bodft_avx2_notes.md

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,3 +37,9 @@
3737
- 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.
3838
- 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.
3939
- 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.
40+
41+
## Float inverse combine pass
42+
43+
- Added an eight-position AVX2/FMA single-precision inverse combine path that mirrors the existing double inverse algebra and reuses the float AVX2 deinterleave/interleave helpers.
44+
- The path loads k/k+M and conjugate-partner blocks, reverses partner lanes in registers, reconstructs the four child spectra, applies conjugate twiddle rotations with `vfmaddps`/`vfmsubps`, and stores c0..c3 as packed complex children.
45+
- In this container, a same-command BODFT benchmark comparison against the previous commit at N=65536 improved float inverse from about 1.13 ms to about 0.78 ms; N=1048576 improved from about 25.6 ms to about 18.4 ms. Forward timings remained noisy and mostly unchanged because this patch targets inverse combine work.

src/detail/bodft_kernel.hpp

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -239,6 +239,70 @@ static inline int bodft_combine_fwd_avx2_f32(const complex_f32_t* RESTRICT tab,
239239
return k;
240240
}
241241

242+
243+
static inline int bodft_combine_inv_avx2_f32(const complex_f32_t* RESTRICT tab,
244+
const complex_f32_t* RESTRICT tab2,
245+
const complex_f32_t* RESTRICT tab3,
246+
const complex_f32_t* RESTRICT in,
247+
complex_f32_t* RESTRICT c0,
248+
complex_f32_t* RESTRICT c1,
249+
complex_f32_t* RESTRICT c2,
250+
complex_f32_t* RESTRICT c3,
251+
int M, int half) {
252+
const __m256 halfv = _mm256_set1_ps(0.5f);
253+
const __m256 zero = _mm256_setzero_ps();
254+
int k = 0;
255+
for (; k + 8 <= half; k += 8) {
256+
__m256 tre, tim, t2re, t2im, t3re, t3im;
257+
__m256 ykre, ykim, ykMre, ykMim;
258+
__m256 ykpre, ykpim, ykpMre, ykpMim;
259+
bodft_uzp8_ps(tab, k, tre, tim);
260+
bodft_uzp8_ps(tab2, k, t2re, t2im);
261+
bodft_uzp8_ps(tab3, k, t3re, t3im);
262+
bodft_uzp8_ps(in, k, ykre, ykim);
263+
bodft_uzp8_ps(in, k + M, ykMre, ykMim);
264+
bodft_uzp8_ps(in, M - 8 - k, ykpre, ykpim);
265+
bodft_uzp8_ps(in, M - 8 - k + M, ykpMre, ykpMim);
266+
ykpre = bodft_rev8_ps(ykpre);
267+
ykpim = bodft_rev8_ps(ykpim);
268+
ykpMre = bodft_rev8_ps(ykpMre);
269+
ykpMim = bodft_rev8_ps(ykpMim);
270+
271+
const __m256 e0re = _mm256_mul_ps(halfv, _mm256_add_ps(ykre, ykpMre));
272+
const __m256 e0im = _mm256_mul_ps(halfv, _mm256_sub_ps(ykim, ykpMim));
273+
const __m256 o0re = _mm256_mul_ps(halfv, _mm256_sub_ps(ykre, ykpMre));
274+
const __m256 o0im = _mm256_mul_ps(halfv, _mm256_add_ps(ykim, ykpMim));
275+
const __m256 e1re = _mm256_mul_ps(halfv, _mm256_add_ps(ykMre, ykpre));
276+
const __m256 e1im = _mm256_mul_ps(halfv, _mm256_sub_ps(ykMim, ykpim));
277+
const __m256 diff_re = _mm256_sub_ps(ykMre, ykpre);
278+
const __m256 diff_im = _mm256_add_ps(ykMim, ykpim);
279+
const __m256 o1re = _mm256_mul_ps(halfv, _mm256_sub_ps(zero, diff_im));
280+
const __m256 o1im = _mm256_mul_ps(halfv, diff_re);
281+
282+
const __m256 b0re = _mm256_mul_ps(halfv, _mm256_add_ps(e0re, e1re));
283+
const __m256 b0im = _mm256_mul_ps(halfv, _mm256_add_ps(e0im, e1im));
284+
const __m256 b2re = _mm256_mul_ps(halfv, _mm256_sub_ps(e0re, e1re));
285+
const __m256 b2im = _mm256_mul_ps(halfv, _mm256_sub_ps(e0im, e1im));
286+
const __m256 b1re = _mm256_mul_ps(halfv, _mm256_add_ps(o0re, o1re));
287+
const __m256 b1im = _mm256_mul_ps(halfv, _mm256_add_ps(o0im, o1im));
288+
const __m256 b3re = _mm256_mul_ps(halfv, _mm256_sub_ps(o0re, o1re));
289+
const __m256 b3im = _mm256_mul_ps(halfv, _mm256_sub_ps(o0im, o1im));
290+
291+
const __m256 c1re = _mm256_fmadd_ps(tim, b1im, _mm256_mul_ps(tre, b1re));
292+
const __m256 c1im = _mm256_fmsub_ps(tre, b1im, _mm256_mul_ps(tim, b1re));
293+
const __m256 c2re = _mm256_fmadd_ps(t2im, b2im, _mm256_mul_ps(t2re, b2re));
294+
const __m256 c2im = _mm256_fmsub_ps(t2re, b2im, _mm256_mul_ps(t2im, b2re));
295+
const __m256 c3re = _mm256_fmadd_ps(t3im, b3im, _mm256_mul_ps(t3re, b3re));
296+
const __m256 c3im = _mm256_fmsub_ps(t3re, b3im, _mm256_mul_ps(t3im, b3re));
297+
298+
bodft_store8_ps(c0, k, b0re, b0im);
299+
bodft_store8_ps(c1, k, c1re, c1im);
300+
bodft_store8_ps(c2, k, c2re, c2im);
301+
bodft_store8_ps(c3, k, c3re, c3im);
302+
}
303+
return k;
304+
}
305+
242306
static inline int bodft_combine_inv_avx2_f64(const complex_t* RESTRICT tab,
243307
const complex_t* RESTRICT tab2,
244308
const complex_t* RESTRICT tab3,
@@ -537,6 +601,17 @@ static inline void combine_inv(const CT* RESTRICT tab, const CT* RESTRICT tab2,
537601
reinterpret_cast<complex_t*>(c2),
538602
reinterpret_cast<complex_t*>(c3), M, half);
539603
}
604+
if constexpr (sizeof(RT) == 4) {
605+
k = bodft_combine_inv_avx2_f32(
606+
reinterpret_cast<const complex_f32_t*>(tab),
607+
reinterpret_cast<const complex_f32_t*>(tab2),
608+
reinterpret_cast<const complex_f32_t*>(tab3),
609+
reinterpret_cast<const complex_f32_t*>(in),
610+
reinterpret_cast<complex_f32_t*>(c0),
611+
reinterpret_cast<complex_f32_t*>(c1),
612+
reinterpret_cast<complex_f32_t*>(c2),
613+
reinterpret_cast<complex_f32_t*>(c3), M, half);
614+
}
540615
#endif
541616
const RT h = static_cast<RT>(0.5);
542617
for (; k < half; ++k) {

src/detail/bruun_kernel.hpp

Lines changed: 70 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -59,8 +59,13 @@ typedef __m128d bruun_v2;
5959
# define V2_SUB(a, b) _mm_sub_pd((a), (b))
6060
# define V2_MUL(a, b) _mm_mul_pd((a), (b))
6161
# define V2_DIV(a, b) _mm_div_pd((a), (b))
62-
# define V2_MADD(a, b, c) V2_ADD((a), V2_MUL((b), (c)))
63-
# define V2_MSUB(a, b, c) V2_SUB((a), V2_MUL((b), (c)))
62+
# if defined(__FMA__)
63+
# define V2_MADD(a, b, c) _mm_fmadd_pd((b), (c), (a))
64+
# define V2_MSUB(a, b, c) _mm_fnmadd_pd((b), (c), (a))
65+
# else
66+
# define V2_MADD(a, b, c) V2_ADD((a), V2_MUL((b), (c)))
67+
# define V2_MSUB(a, b, c) V2_SUB((a), V2_MUL((b), (c)))
68+
# endif
6469
# define V2_SET1(x) _mm_set1_pd(x)
6570
# define V2_SETLH(l, h) _mm_set_pd((h), (l))
6671
# define V2_UNPLO(a, b) _mm_unpacklo_pd((a), (b))
@@ -109,8 +114,13 @@ typedef __m128 bruun_v4f;
109114
# define V4F_ADD(a, b) _mm_add_ps((a), (b))
110115
# define V4F_SUB(a, b) _mm_sub_ps((a), (b))
111116
# define V4F_MUL(a, b) _mm_mul_ps((a), (b))
112-
# define V4F_MADD(a, b, c) V4F_ADD((a), V4F_MUL((b), (c)))
113-
# define V4F_MSUB(a, b, c) V4F_SUB((a), V4F_MUL((b), (c)))
117+
# if defined(__FMA__)
118+
# define V4F_MADD(a, b, c) _mm_fmadd_ps((b), (c), (a))
119+
# define V4F_MSUB(a, b, c) _mm_fnmadd_ps((b), (c), (a))
120+
# else
121+
# define V4F_MADD(a, b, c) V4F_ADD((a), V4F_MUL((b), (c)))
122+
# define V4F_MSUB(a, b, c) V4F_SUB((a), V4F_MUL((b), (c)))
123+
# endif
114124
# define V4F_SET1(x) _mm_set1_ps(x)
115125
# define V4F_SET4(a,b,c,d) _mm_setr_ps((a), (b), (c), (d))
116126
# define V4F_ZERO() _mm_setzero_ps()
@@ -680,7 +690,8 @@ static inline void binomial_fwd(double* RESTRICT v, int h) {
680690
_mm256_storeu_pd(v + i, _mm256_add_pd(a, b));
681691
_mm256_storeu_pd(v + h + i, _mm256_sub_pd(a, b));
682692
}
683-
#elif BRUUN_LEVEL == 1
693+
#endif
694+
#if BRUUN_LEVEL >= 1
684695
for (; i + 1 < h; i += 2) {
685696
const bruun_v2 a = V2_LD(v + i);
686697
const bruun_v2 b = V2_LD(v + h + i);
@@ -716,7 +727,8 @@ static inline void binomial_oop(const double* RESTRICT in, double* RESTRICT v, i
716727
_mm256_storeu_pd(v + i, _mm256_add_pd(a, b));
717728
_mm256_storeu_pd(v + h + i, _mm256_sub_pd(a, b));
718729
}
719-
#elif BRUUN_LEVEL == 1
730+
#endif
731+
#if BRUUN_LEVEL >= 1
720732
for (; i + 1 < h; i += 2) {
721733
const bruun_v2 a = V2_LD(in + i);
722734
const bruun_v2 b = V2_LD(in + h + i);
@@ -755,7 +767,8 @@ static inline void binomial_inv(double* RESTRICT v, int h) {
755767
_mm256_storeu_pd(v + i, _mm256_mul_pd(half4, _mm256_add_pd(a, b)));
756768
_mm256_storeu_pd(v + h + i, _mm256_mul_pd(half4, _mm256_sub_pd(a, b)));
757769
}
758-
#elif BRUUN_LEVEL == 1
770+
#endif
771+
#if BRUUN_LEVEL >= 1
759772
const bruun_v2 half2 = V2_SET1(0.5);
760773
for (; i + 1 < h; i += 2) {
761774
const bruun_v2 a = V2_LD(v + i);
@@ -823,7 +836,8 @@ static inline void norm_q_fwd(double* RESTRICT p, int q, double c_scalar, double
823836
_mm256_storeu_pd(B1p + n, _mm256_sub_pd(I, A1));
824837
}
825838
}
826-
#elif BRUUN_LEVEL == 1
839+
#endif
840+
#if BRUUN_LEVEL >= 1
827841
{
828842
const bruun_v2 vc = V2_SET1(c_scalar);
829843
const bruun_v2 vs = V2_SET1(s_scalar);
@@ -967,7 +981,8 @@ static inline void norm2_fused(double* RESTRICT p, int q,
967981
_mm256_storeu_pd(B1 + qh + n, _mm256_sub_pd(I1, x0));
968982
}
969983
}
970-
#elif BRUUN_LEVEL == 1
984+
#endif
985+
#if BRUUN_LEVEL >= 1
971986
{
972987
const bruun_v2 vc = V2_SET1(c), vs = V2_SET1(s);
973988
const bruun_v2 vc0 = V2_SET1(c0), vs0 = V2_SET1(s0);
@@ -1107,7 +1122,8 @@ static inline void norm_q_inv(double* RESTRICT p, int q, double c_scalar, double
11071122
_mm256_storeu_pd(D1p + n, B1);
11081123
}
11091124
}
1110-
#elif BRUUN_LEVEL == 1
1125+
#endif
1126+
#if BRUUN_LEVEL >= 1
11111127
{
11121128
const bruun_v2 half = V2_SET1(0.5);
11131129
const bruun_v2 vc = V2_SET1(c_scalar);
@@ -1268,7 +1284,8 @@ static inline void norm2_inv_fused(double* RESTRICT p, int q,
12681284
_mm256_storeu_pd(B1 + qh + n, _mm256_fmsub_pd(vc, Ih, _mm256_mul_pd(vs, Rh)));
12691285
}
12701286
}
1271-
#elif BRUUN_LEVEL == 1
1287+
#endif
1288+
#if BRUUN_LEVEL >= 1
12721289
{
12731290
const bruun_v2 hf = V2_SET1(0.5);
12741291
const bruun_v2 vc = V2_SET1(c), vs = V2_SET1(s);
@@ -1532,14 +1549,14 @@ static inline void norm_q_fwd_f32(float* RESTRICT p, int q, float c_scalar, floa
15321549
const __m256 vs = _mm256_set1_ps(s_scalar);
15331550

15341551
for (; n + 7 < q; n += 8) {
1535-
const __m256 A0 = _mm256_loadu_ps(A0p + n);
15361552
const __m256 B0 = _mm256_loadu_ps(B0p + n);
1537-
const __m256 A1 = _mm256_loadu_ps(A1p + n);
15381553
const __m256 B1 = _mm256_loadu_ps(B1p + n);
1539-
15401554
const __m256 R = _mm256_fmsub_ps(vc, B0, _mm256_mul_ps(vs, B1));
15411555
const __m256 I = _mm256_fmadd_ps(vs, B0, _mm256_mul_ps(vc, B1));
15421556

1557+
const __m256 A0 = _mm256_loadu_ps(A0p + n);
1558+
const __m256 A1 = _mm256_loadu_ps(A1p + n);
1559+
15431560
_mm256_storeu_ps(A0p + n, _mm256_add_ps(A0, R));
15441561
_mm256_storeu_ps(B0p + n, _mm256_add_ps(A1, I));
15451562
_mm256_storeu_ps(A1p + n, _mm256_sub_ps(A0, R));
@@ -1553,14 +1570,14 @@ static inline void norm_q_fwd_f32(float* RESTRICT p, int q, float c_scalar, floa
15531570
const bruun_v4f vs = V4F_SET1(s_scalar);
15541571

15551572
for (; n + 3 < q; n += 4) {
1556-
const bruun_v4f A0 = V4F_LD(A0p + n);
15571573
const bruun_v4f B0 = V4F_LD(B0p + n);
1558-
const bruun_v4f A1 = V4F_LD(A1p + n);
15591574
const bruun_v4f B1 = V4F_LD(B1p + n);
1560-
15611575
const bruun_v4f R = V4F_MSUB(V4F_MUL(vc, B0), vs, B1);
15621576
const bruun_v4f I = V4F_MADD(V4F_MUL(vs, B0), vc, B1);
15631577

1578+
const bruun_v4f A0 = V4F_LD(A0p + n);
1579+
const bruun_v4f A1 = V4F_LD(A1p + n);
1580+
15641581
V4F_ST(A0p + n, V4F_ADD(A0, R));
15651582
V4F_ST(B0p + n, V4F_ADD(A1, I));
15661583
V4F_ST(A1p + n, V4F_SUB(A0, R));
@@ -1652,20 +1669,21 @@ static inline void norm2_fused_f32(float* RESTRICT p, int q,
16521669
const __m256 vc1 = _mm256_set1_ps(c1), vs1 = _mm256_set1_ps(s1);
16531670

16541671
for (; n + 7 < qh; n += 8) {
1655-
const __m256 a0n = _mm256_loadu_ps(A0 + n);
1656-
const __m256 a0h = _mm256_loadu_ps(A0 + qh + n);
16571672
const __m256 b0n = _mm256_loadu_ps(B0 + n);
1658-
const __m256 b0h = _mm256_loadu_ps(B0 + qh + n);
1659-
const __m256 a1n = _mm256_loadu_ps(A1 + n);
1660-
const __m256 a1h = _mm256_loadu_ps(A1 + qh + n);
16611673
const __m256 b1n = _mm256_loadu_ps(B1 + n);
1674+
const __m256 b0h = _mm256_loadu_ps(B0 + qh + n);
16621675
const __m256 b1h = _mm256_loadu_ps(B1 + qh + n);
16631676

16641677
const __m256 Rn = _mm256_fmsub_ps(vc, b0n, _mm256_mul_ps(vs, b1n));
16651678
const __m256 In = _mm256_fmadd_ps(vs, b0n, _mm256_mul_ps(vc, b1n));
16661679
const __m256 Rh = _mm256_fmsub_ps(vc, b0h, _mm256_mul_ps(vs, b1h));
16671680
const __m256 Ih = _mm256_fmadd_ps(vs, b0h, _mm256_mul_ps(vc, b1h));
16681681

1682+
const __m256 a0n = _mm256_loadu_ps(A0 + n);
1683+
const __m256 a0h = _mm256_loadu_ps(A0 + qh + n);
1684+
const __m256 a1n = _mm256_loadu_ps(A1 + n);
1685+
const __m256 a1h = _mm256_loadu_ps(A1 + qh + n);
1686+
16691687
const __m256 u0 = _mm256_add_ps(a0n, Rn);
16701688
const __m256 uh = _mm256_add_ps(a0h, Rh);
16711689
const __m256 w0 = _mm256_add_ps(a1n, In);
@@ -1698,20 +1716,21 @@ static inline void norm2_fused_f32(float* RESTRICT p, int q,
16981716
const bruun_v4f vc1 = V4F_SET1(c1), vs1 = V4F_SET1(s1);
16991717

17001718
for (; n + 3 < qh; n += 4) {
1701-
const bruun_v4f a0n = V4F_LD(A0 + n);
1702-
const bruun_v4f a0h = V4F_LD(A0 + qh + n);
17031719
const bruun_v4f b0n = V4F_LD(B0 + n);
1704-
const bruun_v4f b0h = V4F_LD(B0 + qh + n);
1705-
const bruun_v4f a1n = V4F_LD(A1 + n);
1706-
const bruun_v4f a1h = V4F_LD(A1 + qh + n);
17071720
const bruun_v4f b1n = V4F_LD(B1 + n);
1721+
const bruun_v4f b0h = V4F_LD(B0 + qh + n);
17081722
const bruun_v4f b1h = V4F_LD(B1 + qh + n);
17091723

17101724
const bruun_v4f Rn = V4F_MSUB(V4F_MUL(vc, b0n), vs, b1n);
17111725
const bruun_v4f In = V4F_MADD(V4F_MUL(vs, b0n), vc, b1n);
17121726
const bruun_v4f Rh = V4F_MSUB(V4F_MUL(vc, b0h), vs, b1h);
17131727
const bruun_v4f Ih = V4F_MADD(V4F_MUL(vs, b0h), vc, b1h);
17141728

1729+
const bruun_v4f a0n = V4F_LD(A0 + n);
1730+
const bruun_v4f a0h = V4F_LD(A0 + qh + n);
1731+
const bruun_v4f a1n = V4F_LD(A1 + n);
1732+
const bruun_v4f a1h = V4F_LD(A1 + qh + n);
1733+
17151734
const bruun_v4f u0 = V4F_ADD(a0n, Rn);
17161735
const bruun_v4f uh = V4F_ADD(a0h, Rh);
17171736
const bruun_v4f w0 = V4F_ADD(a1n, In);
@@ -1808,22 +1827,24 @@ static inline void norm_q_inv_f32(float* RESTRICT p, int q, float c_scalar, floa
18081827
#if BRUUN_LEVEL >= 2
18091828
{
18101829
const __m256 half = _mm256_set1_ps(0.5f);
1811-
const __m256 vc = _mm256_set1_ps(c_scalar);
1812-
const __m256 vs = _mm256_set1_ps(s_scalar);
1830+
const __m256 hvc = _mm256_set1_ps(0.5f * c_scalar);
1831+
const __m256 hvs = _mm256_set1_ps(0.5f * s_scalar);
18131832

18141833
for (; n + 7 < q; n += 8) {
18151834
const __m256 C0v = _mm256_loadu_ps(C0p + n);
18161835
const __m256 C1v = _mm256_loadu_ps(C1p + n);
18171836
const __m256 D0v = _mm256_loadu_ps(D0p + n);
18181837
const __m256 D1v = _mm256_loadu_ps(D1p + n);
18191838

1820-
const __m256 A0 = _mm256_mul_ps(half, _mm256_add_ps(C0v, D0v));
1821-
const __m256 R = _mm256_mul_ps(half, _mm256_sub_ps(C0v, D0v));
1822-
const __m256 I = _mm256_mul_ps(half, _mm256_add_ps(C1v, D1v));
1823-
const __m256 A1 = _mm256_mul_ps(half, _mm256_sub_ps(C1v, D1v));
1839+
const __m256 t0 = _mm256_add_ps(C0v, D0v);
1840+
const __m256 r = _mm256_sub_ps(C0v, D0v);
1841+
const __m256 i = _mm256_add_ps(C1v, D1v);
1842+
const __m256 t1 = _mm256_sub_ps(C1v, D1v);
18241843

1825-
const __m256 B0 = _mm256_fmadd_ps(vc, R, _mm256_mul_ps(vs, I));
1826-
const __m256 B1 = _mm256_fmsub_ps(vc, I, _mm256_mul_ps(vs, R));
1844+
const __m256 A0 = _mm256_mul_ps(half, t0);
1845+
const __m256 A1 = _mm256_mul_ps(half, t1);
1846+
const __m256 B0 = _mm256_fmadd_ps(hvc, r, _mm256_mul_ps(hvs, i));
1847+
const __m256 B1 = _mm256_fmsub_ps(hvc, i, _mm256_mul_ps(hvs, r));
18271848

18281849
_mm256_storeu_ps(C0p + n, A0);
18291850
_mm256_storeu_ps(C1p + n, B0);
@@ -1835,22 +1856,24 @@ static inline void norm_q_inv_f32(float* RESTRICT p, int q, float c_scalar, floa
18351856
#if BRUUN_LEVEL >= 1
18361857
{
18371858
const bruun_v4f half = V4F_SET1(0.5f);
1838-
const bruun_v4f vc = V4F_SET1(c_scalar);
1839-
const bruun_v4f vs = V4F_SET1(s_scalar);
1859+
const bruun_v4f hvc = V4F_SET1(0.5f * c_scalar);
1860+
const bruun_v4f hvs = V4F_SET1(0.5f * s_scalar);
18401861

18411862
for (; n + 3 < q; n += 4) {
18421863
const bruun_v4f C0v = V4F_LD(C0p + n);
18431864
const bruun_v4f C1v = V4F_LD(C1p + n);
18441865
const bruun_v4f D0v = V4F_LD(D0p + n);
18451866
const bruun_v4f D1v = V4F_LD(D1p + n);
18461867

1847-
const bruun_v4f A0 = V4F_MUL(half, V4F_ADD(C0v, D0v));
1848-
const bruun_v4f R = V4F_MUL(half, V4F_SUB(C0v, D0v));
1849-
const bruun_v4f I = V4F_MUL(half, V4F_ADD(C1v, D1v));
1850-
const bruun_v4f A1 = V4F_MUL(half, V4F_SUB(C1v, D1v));
1868+
const bruun_v4f t0 = V4F_ADD(C0v, D0v);
1869+
const bruun_v4f r = V4F_SUB(C0v, D0v);
1870+
const bruun_v4f i = V4F_ADD(C1v, D1v);
1871+
const bruun_v4f t1 = V4F_SUB(C1v, D1v);
18511872

1852-
const bruun_v4f B0 = V4F_MADD(V4F_MUL(vc, R), vs, I);
1853-
const bruun_v4f B1 = V4F_MSUB(V4F_MUL(vc, I), vs, R);
1873+
const bruun_v4f A0 = V4F_MUL(half, t0);
1874+
const bruun_v4f A1 = V4F_MUL(half, t1);
1875+
const bruun_v4f B0 = V4F_MADD(V4F_MUL(hvc, r), hvs, i);
1876+
const bruun_v4f B1 = V4F_MSUB(V4F_MUL(hvc, i), hvs, r);
18541877

18551878
V4F_ST(C0p + n, A0);
18561879
V4F_ST(C1p + n, B0);
@@ -2367,7 +2390,8 @@ class RFFT {
23672390
bool standard_output_uses_two_phase() const {
23682391
#if BRUUN_LEVEL >= 2
23692392
return N >= 8192;
2370-
#elif BRUUN_LEVEL == 1
2393+
#endif
2394+
#if BRUUN_LEVEL >= 1
23712395
return N > 1048576;
23722396
#else
23732397
return false;
@@ -3045,7 +3069,8 @@ class RFFT {
30453069
_mm256_storeu_pd(v + j, _mm256_fmadd_pd(hd, r, _mm256_mul_pd(ho, rs)));
30463070
}
30473071
}
3048-
#elif BRUUN_LEVEL == 1
3072+
#endif
3073+
#if BRUUN_LEVEL >= 1
30493074
for (; j + 1 < end; j += 2) {
30503075
const bruun_v2 r = V2_LD(v + j);
30513076
const bruun_v2 h = V2_LD(RF + j);

src/detail/inverse_acceleration_notes.txt

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,3 +41,12 @@ Wide-x86 schedule gate:
4141
- Rechecked the scheduled inverse on Intel Xeon Platinum 8370C with `-O3 -mavx2 -mfma -mno-avx512*`.
4242
- Rechecked the AVX-512 scheduled inverse on the same host with `-O3 -march=native`.
4343
- The benchmark prints a direct `Inv_ns` column so inverse-only changes are visible instead of inferred from round-trip timing.
44+
45+
2026-06-27 x86 SIMD tail and f32 scheduling pass:
46+
- x86 128-bit V2/V4F multiply-add macros now use accumulator-first FMA intrinsics when `__FMA__` is available (`_mm_fmadd_*` and `_mm_fnmadd_*`) instead of lowering to separate multiply plus add/subtract. The fallback SSE2 definitions are unchanged for non-FMA builds.
47+
- Double streaming kernels now run their 128-bit V2 loops after the AVX2/AVX-512 loops, so one-to-three-element wide tails no longer drop straight to scalar when a two-lane chunk remains.
48+
- The float forward kernels start the B0/B1 twiddle rotation before loading A0/A1 in the AVX2 and V4F loops, shortening live ranges around the first dependency frontier.
49+
- `norm_q_inv_f32` folds the 0.5 scale into the twiddle constants for rotated B outputs in the AVX2 and V4F loops, keeping explicit half multiplies only for the A outputs that are stored directly.
50+
51+
2026-06-27 continuation:
52+
- Added `kernel_sequentiality_flow.md` to describe forward/inverse dependency frontiers, SIMD tail flow, and where arithmetic can safely move across frontiers without increasing live pressure.

0 commit comments

Comments
 (0)