1#include <assert.h>
 2#include <immintrin.h>
 3#include <stdint.h>
 4
 5#include "reim4_fftvec_private.h"
 6
 7EXPORT void reim4_fftvec_addmul_fma(const REIM4_FFTVEC_ADDMUL_PRECOMP* tables, double* r_ptr, const double* a_ptr,
 8                                    const double* b_ptr) {
 9  const double* const rend_ptr = r_ptr + (tables->m << 1);
10  while (r_ptr != rend_ptr) {
11    __m256d rr = _mm256_loadu_pd(r_ptr);
12    __m256d ri = _mm256_loadu_pd(r_ptr + 4);
13    const __m256d ar = _mm256_loadu_pd(a_ptr);
14    const __m256d ai = _mm256_loadu_pd(a_ptr + 4);
15    const __m256d br = _mm256_loadu_pd(b_ptr);
16    const __m256d bi = _mm256_loadu_pd(b_ptr + 4);
17
18    rr = _mm256_fmsub_pd(ai, bi, rr);
19    rr = _mm256_fmsub_pd(ar, br, rr);
20    ri = _mm256_fmadd_pd(ar, bi, ri);
21    ri = _mm256_fmadd_pd(ai, br, ri);
22
23    _mm256_storeu_pd(r_ptr, rr);
24    _mm256_storeu_pd(r_ptr + 4, ri);
25
26    r_ptr += 8;
27    a_ptr += 8;
28    b_ptr += 8;
29  }
30}
31
32EXPORT void reim4_fftvec_mul_fma(const REIM4_FFTVEC_MUL_PRECOMP* tables, double* r_ptr, const double* a_ptr,
33                                 const double* b_ptr) {
34  const double* const rend_ptr = r_ptr + (tables->m << 1);
35  while (r_ptr != rend_ptr) {
36    const __m256d ar = _mm256_loadu_pd(a_ptr);
37    const __m256d ai = _mm256_loadu_pd(a_ptr + 4);
38    const __m256d br = _mm256_loadu_pd(b_ptr);
39    const __m256d bi = _mm256_loadu_pd(b_ptr + 4);
40
41    const __m256d t1 = _mm256_mul_pd(ai, bi);
42    const __m256d t2 = _mm256_mul_pd(ar, bi);
43
44    __m256d rr = _mm256_fmsub_pd(ar, br, t1);
45    __m256d ri = _mm256_fmadd_pd(ai, br, t2);
46
47    _mm256_storeu_pd(r_ptr, rr);
48    _mm256_storeu_pd(r_ptr + 4, ri);
49
50    r_ptr += 8;
51    a_ptr += 8;
52    b_ptr += 8;
53  }
54}