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(__m256d* ra, __m256d* rb, __m256d* ia, __m256d* ib, const __m256d omre,
 8                                              const __m256d omim) {
 9  __m256d rdiff = _mm256_sub_pd(*ra, *rb);
10  __m256d idiff = _mm256_sub_pd(*ia, *ib);
11  *ra = _mm256_add_pd(*ra, *rb);
12  *ia = _mm256_add_pd(*ia, *ib);
13
14  *rb = _mm256_mul_pd(idiff, omim);
15  *rb = _mm256_fmsub_pd(rdiff, omre, *rb);
16
17  *ib = _mm256_mul_pd(rdiff, omim);
18  *ib = _mm256_fmadd_pd(idiff, omre, *ib);
19}
20
21EXPORT void reim_ifft8_avx_fma(double* dre, double* dim, const void* ompv) {
22  const double* omp = (const double*)ompv;
23
24  __m256d ra0 = _mm256_loadu_pd(dre);
25  __m256d ra4 = _mm256_loadu_pd(dre + 4);
26  __m256d ia0 = _mm256_loadu_pd(dim);
27  __m256d ia4 = _mm256_loadu_pd(dim + 4);
28
29  // 1
30  {
31    const __m256d fft8neg2 = _mm256_castsi256_pd(_mm256_set_epi64x(UINT64_C(1) << 63, UINT64_C(1) << 63, 0, 0));
32    __m256d omr = _mm256_loadu_pd(omp);                    // r0,r1,i0,i1
33    __m256d omiirr = _mm256_permute2f128_pd(omr, omr, 1);  // i0,i1,r0,r1
34    __m256d omi = _mm256_xor_pd(omiirr, fft8neg2);         // i0,i1,-r0,-r1
35
36    __m256d rb = _mm256_unpackhi_pd(ra0, ra4);
37    __m256d ib = _mm256_unpackhi_pd(ia0, ia4);
38    __m256d ra = _mm256_unpacklo_pd(ra0, ra4);
39    __m256d ia = _mm256_unpacklo_pd(ia0, ia4);
40
41    reim_invctwiddle_avx_fma(&ra, &rb, &ia, &ib, omr, omi);
42
43    ra4 = _mm256_unpackhi_pd(ra, rb);
44    ia4 = _mm256_unpackhi_pd(ia, ib);
45    ra0 = _mm256_unpacklo_pd(ra, rb);
46    ia0 = _mm256_unpacklo_pd(ia, ib);
47  }
48
49  // 2
50  {
51    const __m128d ifft8neg = _mm_castsi128_pd(_mm_set_epi64x(0, UINT64_C(1) << 63));
52    __m128d omri = _mm_loadu_pd(omp + 4);             // r,i
53    __m128d ommri = _mm_xor_pd(omri, ifft8neg);       // -r,i
54    __m256d omrimri = _mm256_set_m128d(ommri, omri);  // r,i,-r,i
55    __m256d omi = _mm256_permute_pd(omrimri, 3);      // i,i,-r,-r
56    __m256d omr = _mm256_permute_pd(omrimri, 12);     // r,r,i,i
57
58    __m256d rb = _mm256_permute2f128_pd(ra0, ra4, 0x31);
59    __m256d ib = _mm256_permute2f128_pd(ia0, ia4, 0x31);
60    __m256d ra = _mm256_permute2f128_pd(ra0, ra4, 0x20);
61    __m256d ia = _mm256_permute2f128_pd(ia0, ia4, 0x20);
62
63    reim_invctwiddle_avx_fma(&ra, &rb, &ia, &ib, omr, omi);
64
65    ra0 = _mm256_permute2f128_pd(ra, rb, 0x20);
66    ra4 = _mm256_permute2f128_pd(ra, rb, 0x31);
67    ia0 = _mm256_permute2f128_pd(ia, ib, 0x20);
68    ia4 = _mm256_permute2f128_pd(ia, ib, 0x31);
69  }
70
71  // 3
72  {
73    __m128d omri = _mm_loadu_pd(omp + 6);           // r,i
74    __m256d omriri = _mm256_set_m128d(omri, omri);  // r,i,r,i
75    __m256d omi = _mm256_permute_pd(omriri, 15);    // i,i,i,i
76    __m256d omr = _mm256_permute_pd(omriri, 0);     // r,r,r,r
77
78    reim_invctwiddle_avx_fma(&ra0, &ra4, &ia0, &ia4, omr, omi);
79  }
80
81  // 4
82  _mm256_storeu_pd(dre, ra0);
83  _mm256_storeu_pd(dre + 4, ra4);
84  _mm256_storeu_pd(dim, ia0);
85  _mm256_storeu_pd(dim + 4, ia4);
86}