1#include <inttypes.h>
  2
  3#include <cmath>
  4
  5#include "gtest/gtest.h"
  6#include "spqlios/commons_private.h"
  7#include "spqlios/cplx/cplx_fft_internal.h"
  8#include "spqlios/reim/reim_fft_internal.h"
  9#include "spqlios/reim/reim_fft_private.h"
 10
 11#ifdef __x86_64__
 12TEST(fft, reim_fft_avx2_vs_fft_reim_ref) {
 13  for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
 14    uint64_t m = nn / 2;
 15    // CPLX_FFT_PRECOMP* tables = new_cplx_fft_precomp(m, 0);
 16    REIM_FFT_PRECOMP* reimtables = new_reim_fft_precomp(m, 0);
 17    CPLX* a = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 18    double* a1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 19    double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 20    int64_t p = 1 << 16;
 21    for (uint32_t i = 0; i < nn / 2; i++) {
 22      a[i][0] = (rand() % p) - p / 2;  // between -p/2 and p/2
 23      a[i][1] = (rand() % p) - p / 2;
 24    }
 25    memcpy(a1, a, nn / 2 * sizeof(CPLX));
 26    memcpy(a2, a, nn / 2 * sizeof(CPLX));
 27    reim_fft_ref(reimtables, a2);
 28    reim_fft_avx2_fma(reimtables, a1);
 29    double d = 0;
 30    for (uint32_t i = 0; i < nn / 2; i++) {
 31      double dre = fabs(a1[i] - a2[i]);
 32      double dim = fabs(a1[nn / 2 + i] - a2[nn / 2 + i]);
 33      if (dre > d) d = dre;
 34      if (dim > d) d = dim;
 35      ASSERT_LE(d, nn * 1e-10) << nn;
 36    }
 37    ASSERT_LE(d, nn * 1e-10) << nn;
 38    spqlios_free(a);
 39    spqlios_free(a1);
 40    spqlios_free(a2);
 41    // delete_cplx_fft_precomp(tables);
 42    delete_reim_fft_precomp(reimtables);
 43  }
 44}
 45#endif
 46
 47#ifdef __x86_64__
 48TEST(fft, reim_ifft_avx2_vs_reim_ifft_ref) {
 49  for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
 50    uint64_t m = nn / 2;
 51    // CPLX_FFT_PRECOMP* tables = new_cplx_fft_precomp(m, 0);
 52    REIM_IFFT_PRECOMP* reimtables = new_reim_ifft_precomp(m, 0);
 53    CPLX* a = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 54    double* a1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 55    double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 56    int64_t p = 1 << 16;
 57    for (uint32_t i = 0; i < nn / 2; i++) {
 58      a[i][0] = (rand() % p) - p / 2;  // between -p/2 and p/2
 59      a[i][1] = (rand() % p) - p / 2;
 60    }
 61    memcpy(a1, a, nn / 2 * sizeof(CPLX));
 62    memcpy(a2, a, nn / 2 * sizeof(CPLX));
 63    reim_ifft_ref(reimtables, a2);
 64    reim_ifft_avx2_fma(reimtables, a1);
 65    double d = 0;
 66    for (uint32_t i = 0; i < nn / 2; i++) {
 67      double dre = fabs(a1[i] - a2[i]);
 68      double dim = fabs(a1[nn / 2 + i] - a2[nn / 2 + i]);
 69      if (dre > d) d = dre;
 70      if (dim > d) d = dim;
 71      ASSERT_LE(d, 1e-8);
 72    }
 73    ASSERT_LE(d, 1e-8);
 74    spqlios_free(a);
 75    spqlios_free(a1);
 76    spqlios_free(a2);
 77    // delete_cplx_fft_precomp(tables);
 78    delete_reim_fft_precomp(reimtables);
 79  }
 80}
 81#endif
 82
 83#ifdef __x86_64__
 84TEST(fft, reim_vecfft_addmul_fma_vs_ref) {
 85  for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
 86    uint64_t m = nn / 2;
 87    REIM_FFTVEC_ADDMUL_PRECOMP* tbl = new_reim_fftvec_addmul_precomp(m);
 88    ASSERT_TRUE(tbl != nullptr);
 89    double* a1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 90    double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 91    double* b1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 92    double* b2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 93    double* r1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 94    double* r2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
 95    int64_t p = 1 << 16;
 96    for (uint32_t i = 0; i < nn; i++) {
 97      a1[i] = (rand() % p) - p / 2;  // between -p/2 and p/2
 98      b1[i] = (rand() % p) - p / 2;
 99      r1[i] = (rand() % p) - p / 2;
100    }
101    memcpy(a2, a1, nn / 2 * sizeof(CPLX));
102    memcpy(b2, b1, nn / 2 * sizeof(CPLX));
103    memcpy(r2, r1, nn / 2 * sizeof(CPLX));
104    reim_fftvec_addmul_ref(tbl, r1, a1, b1);
105    reim_fftvec_addmul_fma(tbl, r2, a2, b2);
106    double d = 0;
107    for (uint32_t i = 0; i < nn; i++) {
108      double di = fabs(r1[i] - r2[i]);
109      if (di > d) d = di;
110      ASSERT_LE(d, 1e-8);
111    }
112    ASSERT_LE(d, 1e-8);
113    spqlios_free(a1);
114    spqlios_free(a2);
115    spqlios_free(b1);
116    spqlios_free(b2);
117    spqlios_free(r1);
118    spqlios_free(r2);
119    delete_reim_fftvec_addmul_precomp(tbl);
120  }
121}
122#endif
123
124#ifdef __x86_64__
125TEST(fft, reim_vecfft_mul_fma_vs_ref) {
126  for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
127    uint64_t m = nn / 2;
128    REIM_FFTVEC_MUL_PRECOMP* tbl = new_reim_fftvec_mul_precomp(m);
129    double* a1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
130    double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
131    double* b1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
132    double* b2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
133    double* r1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
134    double* r2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
135    int64_t p = 1 << 16;
136    for (uint32_t i = 0; i < nn; i++) {
137      a1[i] = (rand() % p) - p / 2;  // between -p/2 and p/2
138      b1[i] = (rand() % p) - p / 2;
139      r1[i] = (rand() % p) - p / 2;
140    }
141    memcpy(a2, a1, nn / 2 * sizeof(CPLX));
142    memcpy(b2, b1, nn / 2 * sizeof(CPLX));
143    memcpy(r2, r1, nn / 2 * sizeof(CPLX));
144    reim_fftvec_mul_ref(tbl, r1, a1, b1);
145    reim_fftvec_mul_fma(tbl, r2, a2, b2);
146    double d = 0;
147    for (uint32_t i = 0; i < nn; i++) {
148      double di = fabs(r1[i] - r2[i]);
149      if (di > d) d = di;
150      ASSERT_LE(d, 1e-8);
151    }
152    ASSERT_LE(d, 1e-8);
153    spqlios_free(a1);
154    spqlios_free(a2);
155    spqlios_free(b1);
156    spqlios_free(b2);
157    spqlios_free(r1);
158    spqlios_free(r2);
159    delete_reim_fftvec_mul_precomp(tbl);
160  }
161}
162#endif
163
164#ifdef __x86_64__
165TEST(fft, reim_vecfft_add_fma_vs_ref) {
166  for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
167    uint64_t m = nn / 2;
168    REIM_FFTVEC_ADD_PRECOMP* tbl = new_reim_fftvec_add_precomp(m);
169    double* a1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
170    double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
171    double* b1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
172    double* b2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
173    double* r1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
174    double* r2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
175    int64_t p = 1 << 16;
176    for (uint32_t i = 0; i < nn; i++) {
177      a1[i] = (rand() % p) - p / 2;  // between -p/2 and p/2
178      b1[i] = (rand() % p) - p / 2;
179      r1[i] = (rand() % p) - p / 2;
180    }
181    memcpy(a2, a1, nn / 2 * sizeof(CPLX));
182    memcpy(b2, b1, nn / 2 * sizeof(CPLX));
183    memcpy(r2, r1, nn / 2 * sizeof(CPLX));
184    reim_fftvec_add_ref(tbl, r1, a1, b1);
185    reim_fftvec_add_fma(tbl, r2, a2, b2);
186    double d = 0;
187    for (uint32_t i = 0; i < nn; i++) {
188      double di = fabs(r1[i] - r2[i]);
189      if (di > d) d = di;
190      ASSERT_LE(d, 1e-8);
191    }
192    ASSERT_LE(d, 1e-8);
193    spqlios_free(a1);
194    spqlios_free(a2);
195    spqlios_free(b1);
196    spqlios_free(b2);
197    spqlios_free(r1);
198    spqlios_free(r2);
199    delete_reim_fftvec_add_precomp(tbl);
200  }
201}
202#endif
203
204typedef void (*FILL_REIM_FFT_OMG_F)(const double entry_pwr, double** omg);
205typedef void (*REIM_FFT_F)(double* dre, double* dim, const void* omega);
206
207// template to test a fixed-dimension fft vs. naive
208template <uint64_t N>
209void test_reim_fft_ref_vs_naive(FILL_REIM_FFT_OMG_F fill_omega_f, REIM_FFT_F reim_fft_f) {
210  double om[N];
211  double data[2 * N];
212  double datacopy[2 * N];
213  double* omg = om;
214  fill_omega_f(0.25, &omg);
215  ASSERT_EQ(omg - om, ptrdiff_t(N));  // it may depend on N
216  for (uint64_t i = 0; i < N; ++i) {
217    datacopy[i] = data[i] = (rand() % 100) - 50;
218    datacopy[N + i] = data[N + i] = (rand() % 100) - 50;
219  }
220  reim_fft_f(datacopy, datacopy + N, om);
221  reim_naive_fft(N, 0.25, data, data + N);
222  double d = 0;
223  for (uint64_t i = 0; i < 2 * N; ++i) {
224    d += fabs(datacopy[i] - data[i]);
225  }
226  ASSERT_LE(d, 1e-7);
227}
228
229template <uint64_t N>
230void test_reim_fft_ref_vs_accel(REIM_FFT_F reim_fft_ref_f, REIM_FFT_F reim_fft_accel_f) {
231  double om[N];
232  double data[2 * N];
233  double datacopy[2 * N];
234  for (uint64_t i = 0; i < N; ++i) {
235    om[i] = (rand() % 100) - 50;
236    datacopy[i] = data[i] = (rand() % 100) - 50;
237    datacopy[N + i] = data[N + i] = (rand() % 100) - 50;
238  }
239  reim_fft_ref_f(datacopy, datacopy + N, om);
240  reim_fft_accel_f(data, data + N, om);
241  double d = 0;
242  for (uint64_t i = 0; i < 2 * N; ++i) {
243    d += fabs(datacopy[i] - data[i]);
244  }
245  if (d > 1e-15) {
246    for (uint64_t i = 0; i < N; ++i) {
247      printf("%" PRId64 " %lf %lf %lf %lf\n", i, data[i], data[N + i], datacopy[i], datacopy[N + i]);
248    }
249    ASSERT_LE(d, 0);
250  }
251}
252
253TEST(fft, reim_fft16_ref_vs_naive) { test_reim_fft_ref_vs_naive<16>(fill_reim_fft16_omegas, reim_fft16_ref); }
254#ifdef __aarch64__
255TEST(fft, reim_fft16_neon_vs_naive) { test_reim_fft_ref_vs_naive<16>(fill_reim_fft16_omegas_neon, reim_fft16_neon); }
256#endif
257
258#ifdef __x86_64__
259TEST(fft, reim_fft16_ref_vs_fma) { test_reim_fft_ref_vs_accel<16>(reim_fft16_ref, reim_fft16_avx_fma); }
260#endif
261
262#ifdef __aarch64__
263static void reim_fft16_ref_neon_pom(double* dre, double* dim, const void* omega) {
264  const double* pom = (double*)omega;
265  // put the omegas in neon order
266  double x_pom[] = {pom[0], pom[1],  pom[2],  pom[3],  pom[4], pom[5],  pom[6],  pom[7],
267                    pom[8], pom[10], pom[12], pom[14], pom[9], pom[11], pom[13], pom[15]};
268  reim_fft16_ref(dre, dim, x_pom);
269}
270TEST(fft, reim_fft16_ref_vs_neon) { test_reim_fft_ref_vs_accel<16>(reim_fft16_ref_neon_pom, reim_fft16_neon); }
271#endif
272
273TEST(fft, reim_fft8_ref_vs_naive) { test_reim_fft_ref_vs_naive<8>(fill_reim_fft8_omegas, reim_fft8_ref); }
274
275#ifdef __x86_64__
276TEST(fft, reim_fft8_ref_vs_fma) { test_reim_fft_ref_vs_accel<8>(reim_fft8_ref, reim_fft8_avx_fma); }
277#endif
278
279TEST(fft, reim_fft4_ref_vs_naive) { test_reim_fft_ref_vs_naive<4>(fill_reim_fft4_omegas, reim_fft4_ref); }
280
281#ifdef __x86_64__
282TEST(fft, reim_fft4_ref_vs_fma) { test_reim_fft_ref_vs_accel<4>(reim_fft4_ref, reim_fft4_avx_fma); }
283#endif
284
285TEST(fft, reim_fft2_ref_vs_naive) { test_reim_fft_ref_vs_naive<2>(fill_reim_fft2_omegas, reim_fft2_ref); }
286
287TEST(fft, reim_fft_bfs_16_ref_vs_naive) {
288  for (const uint64_t m : {16, 32, 64, 128, 256, 512, 1024, 2048}) {
289    std::vector<double> om(2 * m);
290    std::vector<double> data(2 * m);
291    std::vector<double> datacopy(2 * m);
292    double* omg = om.data();
293    fill_reim_fft_bfs_16_omegas(m, 0.25, &omg);
294    ASSERT_LE(omg - om.data(), ptrdiff_t(2 * m));  // it may depend on m
295    for (uint64_t i = 0; i < m; ++i) {
296      datacopy[i] = data[i] = (rand() % 100) - 50;
297      datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
298    }
299    omg = om.data();
300    reim_fft_bfs_16_ref(m, datacopy.data(), datacopy.data() + m, &omg);
301    reim_naive_fft(m, 0.25, data.data(), data.data() + m);
302    double d = 0;
303    for (uint64_t i = 0; i < 2 * m; ++i) {
304      d += fabs(datacopy[i] - data[i]);
305    }
306    ASSERT_LE(d, 1e-7);
307  }
308}
309
310TEST(fft, reim_fft_rec_16_ref_vs_naive) {
311  for (const uint64_t m : {2048, 4096, 8192, 32768, 65536}) {
312    std::vector<double> om(2 * m);
313    std::vector<double> data(2 * m);
314    std::vector<double> datacopy(2 * m);
315    double* omg = om.data();
316    fill_reim_fft_rec_16_omegas(m, 0.25, &omg);
317    ASSERT_LE(omg - om.data(), ptrdiff_t(2 * m));  // it may depend on m
318    for (uint64_t i = 0; i < m; ++i) {
319      datacopy[i] = data[i] = (rand() % 100) - 50;
320      datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
321    }
322    omg = om.data();
323    reim_fft_rec_16_ref(m, datacopy.data(), datacopy.data() + m, &omg);
324    reim_naive_fft(m, 0.25, data.data(), data.data() + m);
325    double d = 0;
326    for (uint64_t i = 0; i < 2 * m; ++i) {
327      d += fabs(datacopy[i] - data[i]);
328    }
329    ASSERT_LE(d, 1e-5);
330  }
331}
332
333TEST(fft, reim_fft_ref_vs_naive) {
334  for (const uint64_t m : {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 32768, 65536}) {
335    std::vector<double> om(2 * m);
336    std::vector<double> data(2 * m);
337    std::vector<double> datacopy(2 * m);
338    REIM_FFT_PRECOMP* precomp = new_reim_fft_precomp(m, 0);
339    for (uint64_t i = 0; i < m; ++i) {
340      datacopy[i] = data[i] = (rand() % 100) - 50;
341      datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
342    }
343    reim_fft_ref(precomp, datacopy.data());
344    reim_naive_fft(m, 0.25, data.data(), data.data() + m);
345    double d = 0;
346    for (uint64_t i = 0; i < 2 * m; ++i) {
347      d += fabs(datacopy[i] - data[i]);
348    }
349    ASSERT_LE(d, 1e-5) << m;
350    delete_reim_fft_precomp(precomp);
351  }
352}
353
354#ifdef __aarch64__
355EXPORT REIM_FFT_PRECOMP* new_reim_fft_precomp_neon(uint32_t m, uint32_t num_buffers);
356EXPORT void reim_fft_neon(const REIM_FFT_PRECOMP* precomp, double* d);
357TEST(fft, reim_fft_neon_vs_naive) {
358  for (const uint64_t m : {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 32768, 65536}) {
359    std::vector<double> om(2 * m);
360    std::vector<double> data(2 * m);
361    std::vector<double> datacopy(2 * m);
362    REIM_FFT_PRECOMP* precomp = new_reim_fft_precomp_neon(m, 0);
363    for (uint64_t i = 0; i < m; ++i) {
364      datacopy[i] = data[i] = (rand() % 100) - 50;
365      datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
366    }
367    reim_fft_neon(precomp, datacopy.data());
368    reim_naive_fft(m, 0.25, data.data(), data.data() + m);
369    double d = 0;
370    for (uint64_t i = 0; i < 2 * m; ++i) {
371      d += fabs(datacopy[i] - data[i]);
372    }
373    ASSERT_LE(d, 1e-5) << m;
374    delete_reim_fft_precomp(precomp);
375  }
376}
377#endif
378
379typedef void (*FILL_REIM_IFFT_OMG_F)(const double entry_pwr, double** omg);
380typedef void (*REIM_IFFT_F)(double* dre, double* dim, const void* omega);
381
382// template to test a fixed-dimension fft vs. naive
383template <uint64_t N>
384void test_reim_ifft_ref_vs_naive(FILL_REIM_IFFT_OMG_F fill_omega_f, REIM_IFFT_F reim_ifft_f) {
385  double om[N];
386  double data[2 * N];
387  double datacopy[2 * N];
388  double* omg = om;
389  fill_omega_f(0.25, &omg);
390  ASSERT_EQ(omg - om, ptrdiff_t(N));  // it may depend on N
391  for (uint64_t i = 0; i < N; ++i) {
392    datacopy[i] = data[i] = (rand() % 100) - 50;
393    datacopy[N + i] = data[N + i] = (rand() % 100) - 50;
394  }
395  reim_ifft_f(datacopy, datacopy + N, om);
396  reim_naive_ifft(N, 0.25, data, data + N);
397  double d = 0;
398  for (uint64_t i = 0; i < 2 * N; ++i) {
399    d += fabs(datacopy[i] - data[i]);
400  }
401  ASSERT_LE(d, 1e-7);
402}
403
404template <uint64_t N>
405void test_reim_ifft_ref_vs_accel(REIM_IFFT_F reim_ifft_ref_f, REIM_IFFT_F reim_ifft_accel_f) {
406  double om[N];
407  double data[2 * N];
408  double datacopy[2 * N];
409  for (uint64_t i = 0; i < N; ++i) {
410    om[i] = (rand() % 100) - 50;
411    datacopy[i] = data[i] = (rand() % 100) - 50;
412    datacopy[N + i] = data[N + i] = (rand() % 100) - 50;
413  }
414  reim_ifft_ref_f(datacopy, datacopy + N, om);
415  reim_ifft_accel_f(data, data + N, om);
416  double d = 0;
417  for (uint64_t i = 0; i < 2 * N; ++i) {
418    d += fabs(datacopy[i] - data[i]);
419  }
420  if (d > 1e-15) {
421    for (uint64_t i = 0; i < N; ++i) {
422      printf("%" PRId64 " %lf %lf %lf %lf\n", i, data[i], data[N + i], datacopy[i], datacopy[N + i]);
423    }
424    ASSERT_LE(d, 0);
425  }
426}
427
428TEST(fft, reim_ifft16_ref_vs_naive) { test_reim_ifft_ref_vs_naive<16>(fill_reim_ifft16_omegas, reim_ifft16_ref); }
429
430#ifdef __x86_64__
431TEST(fft, reim_ifft16_ref_vs_fma) { test_reim_ifft_ref_vs_accel<16>(reim_ifft16_ref, reim_ifft16_avx_fma); }
432#endif
433
434TEST(fft, reim_ifft8_ref_vs_naive) { test_reim_ifft_ref_vs_naive<8>(fill_reim_ifft8_omegas, reim_ifft8_ref); }
435
436#ifdef __x86_64__
437TEST(fft, reim_ifft8_ref_vs_fma) { test_reim_ifft_ref_vs_accel<8>(reim_ifft8_ref, reim_ifft8_avx_fma); }
438#endif
439
440TEST(fft, reim_ifft4_ref_vs_naive) { test_reim_ifft_ref_vs_naive<4>(fill_reim_ifft4_omegas, reim_ifft4_ref); }
441
442#ifdef __x86_64__
443TEST(fft, reim_ifft4_ref_vs_fma) { test_reim_ifft_ref_vs_accel<4>(reim_ifft4_ref, reim_ifft4_avx_fma); }
444#endif
445
446TEST(fft, reim_ifft2_ref_vs_naive) { test_reim_ifft_ref_vs_naive<2>(fill_reim_ifft2_omegas, reim_ifft2_ref); }
447
448TEST(fft, reim_ifft_bfs_16_ref_vs_naive) {
449  for (const uint64_t m : {16, 32, 64, 128, 256, 512, 1024, 2048}) {
450    std::vector<double> om(2 * m);
451    std::vector<double> data(2 * m);
452    std::vector<double> datacopy(2 * m);
453    double* omg = om.data();
454    fill_reim_ifft_bfs_16_omegas(m, 0.25, &omg);
455    ASSERT_LE(omg - om.data(), ptrdiff_t(2 * m));  // it may depend on m
456    for (uint64_t i = 0; i < m; ++i) {
457      datacopy[i] = data[i] = (rand() % 100) - 50;
458      datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
459    }
460    omg = om.data();
461    reim_ifft_bfs_16_ref(m, datacopy.data(), datacopy.data() + m, &omg);
462    reim_naive_ifft(m, 0.25, data.data(), data.data() + m);
463    double d = 0;
464    for (uint64_t i = 0; i < 2 * m; ++i) {
465      d += fabs(datacopy[i] - data[i]);
466    }
467    ASSERT_LE(d, 1e-7);
468  }
469}
470
471TEST(fft, reim_ifft_rec_16_ref_vs_naive) {
472  for (const uint64_t m : {2048, 4096, 8192, 32768, 65536}) {
473    std::vector<double> om(2 * m);
474    std::vector<double> data(2 * m);
475    std::vector<double> datacopy(2 * m);
476    double* omg = om.data();
477    fill_reim_ifft_rec_16_omegas(m, 0.25, &omg);
478    ASSERT_LE(omg - om.data(), ptrdiff_t(2 * m));  // it may depend on m
479    for (uint64_t i = 0; i < m; ++i) {
480      datacopy[i] = data[i] = (rand() % 100) - 50;
481      datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
482    }
483    omg = om.data();
484    reim_ifft_rec_16_ref(m, datacopy.data(), datacopy.data() + m, &omg);
485    reim_naive_ifft(m, 0.25, data.data(), data.data() + m);
486    double d = 0;
487    for (uint64_t i = 0; i < 2 * m; ++i) {
488      d += fabs(datacopy[i] - data[i]);
489    }
490    ASSERT_LE(d, 1e-5);
491  }
492}
493
494TEST(fft, reim_ifft_ref_vs_naive) {
495  for (const uint64_t m : {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 32768, 65536}) {
496    std::vector<double> om(2 * m);
497    std::vector<double> data(2 * m);
498    std::vector<double> datacopy(2 * m);
499    REIM_IFFT_PRECOMP* precomp = new_reim_ifft_precomp(m, 0);
500    for (uint64_t i = 0; i < m; ++i) {
501      datacopy[i] = data[i] = (rand() % 100) - 50;
502      datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
503    }
504    reim_ifft_ref(precomp, datacopy.data());
505    reim_naive_ifft(m, 0.25, data.data(), data.data() + m);
506    double d = 0;
507    for (uint64_t i = 0; i < 2 * m; ++i) {
508      d += fabs(datacopy[i] - data[i]);
509    }
510    ASSERT_LE(d, 1e-5) << m;
511    delete_reim_ifft_precomp(precomp);
512  }
513}