1#include <assert.h>
 2#include <immintrin.h>
 3#include <stdint.h>
 4
 5#include "reim_fft_private.h"
 6
 7__always_inline void reim_invctwiddle_avx_fma(__m128d* ra, __m128d* rb, __m128d* ia, __m128d* ib, const __m128d omre,
 8                                              const __m128d omim) {
 9  __m128d rdiff = _mm_sub_pd(*ra, *rb);
10  __m128d idiff = _mm_sub_pd(*ia, *ib);
11  *ra = _mm_add_pd(*ra, *rb);
12  *ia = _mm_add_pd(*ia, *ib);
13
14  *rb = _mm_mul_pd(idiff, omim);
15  *rb = _mm_fmsub_pd(rdiff, omre, *rb);
16
17  *ib = _mm_mul_pd(rdiff, omim);
18  *ib = _mm_fmadd_pd(idiff, omre, *ib);
19}
20
21EXPORT void reim_ifft4_avx_fma(double* dre, double* dim, const void* ompv) {
22  const double* omp = (const double*)ompv;
23
24  __m128d ra0 = _mm_loadu_pd(dre);
25  __m128d ra2 = _mm_loadu_pd(dre + 2);
26  __m128d ia0 = _mm_loadu_pd(dim);
27  __m128d ia2 = _mm_loadu_pd(dim + 2);
28
29  // 1
30  {
31    const __m128d ifft4neg = _mm_castsi128_pd(_mm_set_epi64x(UINT64_C(1) << 63, 0));
32    __m128d omre = _mm_loadu_pd(omp);                              // omre: r,i
33    __m128d omim = _mm_xor_pd(_mm_permute_pd(omre, 1), ifft4neg);  // omim: i,-r
34
35    __m128d ra = _mm_unpacklo_pd(ra0, ra2);  // (r0, r1), (r2, r3) -> (r0, r2)
36    __m128d ia = _mm_unpacklo_pd(ia0, ia2);  // (i0, i1), (i2, i3) -> (i0, i2)
37    __m128d rb = _mm_unpackhi_pd(ra0, ra2);  // (r0, r1), (r2, r3) -> (r1, r3)
38    __m128d ib = _mm_unpackhi_pd(ia0, ia2);  // (i0, i1), (i2, i3) -> (i1, i3)
39
40    reim_invctwiddle_avx_fma(&ra, &rb, &ia, &ib, omre, omim);
41
42    ra0 = _mm_unpacklo_pd(ra, rb);
43    ia0 = _mm_unpacklo_pd(ia, ib);
44    ra2 = _mm_unpackhi_pd(ra, rb);
45    ia2 = _mm_unpackhi_pd(ia, ib);
46  }
47
48  // 2
49  {
50    __m128d om = _mm_loadu_pd(omp + 2);
51    __m128d omre = _mm_permute_pd(om, 0);
52    __m128d omim = _mm_permute_pd(om, 3);
53
54    reim_invctwiddle_avx_fma(&ra0, &ra2, &ia0, &ia2, omre, omim);
55  }
56
57  // 4
58  _mm_storeu_pd(dre, ra0);
59  _mm_storeu_pd(dre + 2, ra2);
60  _mm_storeu_pd(dim, ia0);
61  _mm_storeu_pd(dim + 2, ia2);
62}