1#include <assert.h>
  2#include <immintrin.h>
  3#include <stdint.h>
  4
  5#include "reim_fft_private.h"
  6
  7EXPORT void reim_fftvec_addmul_fma(const REIM_FFTVEC_ADDMUL_PRECOMP* precomp, double* r, const double* a,
  8                                   const double* b) {
  9  const uint64_t m = precomp->m;
 10  double* rr_ptr = r;
 11  double* ri_ptr = r + m;
 12  const double* ar_ptr = a;
 13  const double* ai_ptr = a + m;
 14  const double* br_ptr = b;
 15  const double* bi_ptr = b + m;
 16
 17  const double* const rend_ptr = ri_ptr;
 18  while (rr_ptr != rend_ptr) {
 19    __m256d rr = _mm256_loadu_pd(rr_ptr);
 20    __m256d ri = _mm256_loadu_pd(ri_ptr);
 21    const __m256d ar = _mm256_loadu_pd(ar_ptr);
 22    const __m256d ai = _mm256_loadu_pd(ai_ptr);
 23    const __m256d br = _mm256_loadu_pd(br_ptr);
 24    const __m256d bi = _mm256_loadu_pd(bi_ptr);
 25
 26    rr = _mm256_fmsub_pd(ai, bi, rr);
 27    rr = _mm256_fmsub_pd(ar, br, rr);
 28    ri = _mm256_fmadd_pd(ar, bi, ri);
 29    ri = _mm256_fmadd_pd(ai, br, ri);
 30
 31    _mm256_storeu_pd(rr_ptr, rr);
 32    _mm256_storeu_pd(ri_ptr, ri);
 33
 34    rr_ptr += 4;
 35    ri_ptr += 4;
 36    ar_ptr += 4;
 37    ai_ptr += 4;
 38    br_ptr += 4;
 39    bi_ptr += 4;
 40  }
 41}
 42
 43EXPORT void reim_fftvec_mul_fma(const REIM_FFTVEC_MUL_PRECOMP* precomp, double* r, const double* a, const double* b) {
 44  const uint64_t m = precomp->m;
 45  double* rr_ptr = r;
 46  double* ri_ptr = r + m;
 47  const double* ar_ptr = a;
 48  const double* ai_ptr = a + m;
 49  const double* br_ptr = b;
 50  const double* bi_ptr = b + m;
 51
 52  const double* const rend_ptr = ri_ptr;
 53  while (rr_ptr != rend_ptr) {
 54    const __m256d ar = _mm256_loadu_pd(ar_ptr);
 55    const __m256d ai = _mm256_loadu_pd(ai_ptr);
 56    const __m256d br = _mm256_loadu_pd(br_ptr);
 57    const __m256d bi = _mm256_loadu_pd(bi_ptr);
 58
 59    const __m256d t1 = _mm256_mul_pd(ai, bi);
 60    const __m256d t2 = _mm256_mul_pd(ar, bi);
 61
 62    __m256d rr = _mm256_fmsub_pd(ar, br, t1);
 63    __m256d ri = _mm256_fmadd_pd(ai, br, t2);
 64
 65    _mm256_storeu_pd(rr_ptr, rr);
 66    _mm256_storeu_pd(ri_ptr, ri);
 67
 68    rr_ptr += 4;
 69    ri_ptr += 4;
 70    ar_ptr += 4;
 71    ai_ptr += 4;
 72    br_ptr += 4;
 73    bi_ptr += 4;
 74  }
 75}
 76
 77EXPORT void reim_fftvec_add_fma(const REIM_FFTVEC_ADD_PRECOMP* precomp, double* r, const double* a, const double* b) {
 78  const uint64_t m = precomp->m;
 79  double* rr_ptr = r;
 80  double* ri_ptr = r + m;
 81  const double* ar_ptr = a;
 82  const double* ai_ptr = a + m;
 83  const double* br_ptr = b;
 84  const double* bi_ptr = b + m;
 85
 86  const double* const rend_ptr = ri_ptr;
 87  while (rr_ptr != rend_ptr) {
 88    const __m256d ar = _mm256_loadu_pd(ar_ptr);
 89    const __m256d ai = _mm256_loadu_pd(ai_ptr);
 90    const __m256d br = _mm256_loadu_pd(br_ptr);
 91    const __m256d bi = _mm256_loadu_pd(bi_ptr);
 92
 93    __m256d rr = _mm256_add_pd(ar, br);
 94    __m256d ri = _mm256_add_pd(ai, bi);
 95
 96    _mm256_storeu_pd(rr_ptr, rr);
 97    _mm256_storeu_pd(ri_ptr, ri);
 98
 99    rr_ptr += 4;
100    ri_ptr += 4;
101    ar_ptr += 4;
102    ai_ptr += 4;
103    br_ptr += 4;
104    bi_ptr += 4;
105  }
106}
107
108EXPORT void reim_fftvec_sub_fma(const REIM_FFTVEC_ADD_PRECOMP* precomp, double* r, const double* a, const double* b) {
109  const uint64_t m = precomp->m;
110  double* rr_ptr = r;
111  double* ri_ptr = r + m;
112  const double* ar_ptr = a;
113  const double* ai_ptr = a + m;
114  const double* br_ptr = b;
115  const double* bi_ptr = b + m;
116
117  const double* const rend_ptr = ri_ptr;
118  while (rr_ptr != rend_ptr) {
119    const __m256d ar = _mm256_loadu_pd(ar_ptr);
120    const __m256d ai = _mm256_loadu_pd(ai_ptr);
121    const __m256d br = _mm256_loadu_pd(br_ptr);
122    const __m256d bi = _mm256_loadu_pd(bi_ptr);
123
124    __m256d rr = _mm256_sub_pd(ar, br);
125    __m256d ri = _mm256_sub_pd(ai, bi);
126
127    _mm256_storeu_pd(rr_ptr, rr);
128    _mm256_storeu_pd(ri_ptr, ri);
129
130    rr_ptr += 4;
131    ri_ptr += 4;
132    ar_ptr += 4;
133    ai_ptr += 4;
134    br_ptr += 4;
135    bi_ptr += 4;
136  }
137}