poulpy/spqlios/lib/test/spqlios_reim_test.cpp
spqlios/lib/test/spqlios_reim_test.cppRaw
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
164typedef void (*FILL_REIM_FFT_OMG_F)(const double entry_pwr, double** omg);
165typedef void (*REIM_FFT_F)(double* dre, double* dim, const void* omega);
166
167// template to test a fixed-dimension fft vs. naive
168template <uint64_t N>
169void test_reim_fft_ref_vs_naive(FILL_REIM_FFT_OMG_F fill_omega_f, REIM_FFT_F reim_fft_f) {
170 double om[N];
171 double data[2 * N];
172 double datacopy[2 * N];
173 double* omg = om;
174 fill_omega_f(0.25, &omg);
175 ASSERT_EQ(omg - om, ptrdiff_t(N)); // it may depend on N
176 for (uint64_t i = 0; i < N; ++i) {
177 datacopy[i] = data[i] = (rand() % 100) - 50;
178 datacopy[N + i] = data[N + i] = (rand() % 100) - 50;
179 }
180 reim_fft_f(datacopy, datacopy + N, om);
181 reim_naive_fft(N, 0.25, data, data + N);
182 double d = 0;
183 for (uint64_t i = 0; i < 2 * N; ++i) {
184 d += fabs(datacopy[i] - data[i]);
185 }
186 ASSERT_LE(d, 1e-7);
187}
188
189template <uint64_t N>
190void test_reim_fft_ref_vs_accel(REIM_FFT_F reim_fft_ref_f, REIM_FFT_F reim_fft_accel_f) {
191 double om[N];
192 double data[2 * N];
193 double datacopy[2 * N];
194 for (uint64_t i = 0; i < N; ++i) {
195 om[i] = (rand() % 100) - 50;
196 datacopy[i] = data[i] = (rand() % 100) - 50;
197 datacopy[N + i] = data[N + i] = (rand() % 100) - 50;
198 }
199 reim_fft_ref_f(datacopy, datacopy + N, om);
200 reim_fft_accel_f(data, data + N, om);
201 double d = 0;
202 for (uint64_t i = 0; i < 2 * N; ++i) {
203 d += fabs(datacopy[i] - data[i]);
204 }
205 if (d > 1e-15) {
206 for (uint64_t i = 0; i < N; ++i) {
207 printf("%" PRId64 " %lf %lf %lf %lf\n", i, data[i], data[N + i], datacopy[i], datacopy[N + i]);
208 }
209 ASSERT_LE(d, 0);
210 }
211}
212
213TEST(fft, reim_fft16_ref_vs_naive) { test_reim_fft_ref_vs_naive<16>(fill_reim_fft16_omegas, reim_fft16_ref); }
214#ifdef __aarch64__
215TEST(fft, reim_fft16_neon_vs_naive) { test_reim_fft_ref_vs_naive<16>(fill_reim_fft16_omegas_neon, reim_fft16_neon); }
216#endif
217
218#ifdef __x86_64__
219TEST(fft, reim_fft16_ref_vs_fma) { test_reim_fft_ref_vs_accel<16>(reim_fft16_ref, reim_fft16_avx_fma); }
220#endif
221
222#ifdef __aarch64__
223static void reim_fft16_ref_neon_pom(double* dre, double* dim, const void* omega) {
224 const double* pom = (double*) omega;
225 // put the omegas in neon order
226 double x_pom[] = {
227 pom[0], pom[1], pom[2], pom[3],
228 pom[4],pom[5], pom[6], pom[7],
229 pom[8], pom[10],pom[12], pom[14],
230 pom[9], pom[11],pom[13], pom[15]
231 };
232 reim_fft16_ref(dre, dim, x_pom);
233}
234TEST(fft, reim_fft16_ref_vs_neon) { test_reim_fft_ref_vs_accel<16>(reim_fft16_ref_neon_pom, reim_fft16_neon); }
235#endif
236
237TEST(fft, reim_fft8_ref_vs_naive) { test_reim_fft_ref_vs_naive<8>(fill_reim_fft8_omegas, reim_fft8_ref); }
238
239#ifdef __x86_64__
240TEST(fft, reim_fft8_ref_vs_fma) { test_reim_fft_ref_vs_accel<8>(reim_fft8_ref, reim_fft8_avx_fma); }
241#endif
242
243TEST(fft, reim_fft4_ref_vs_naive) { test_reim_fft_ref_vs_naive<4>(fill_reim_fft4_omegas, reim_fft4_ref); }
244
245#ifdef __x86_64__
246TEST(fft, reim_fft4_ref_vs_fma) { test_reim_fft_ref_vs_accel<4>(reim_fft4_ref, reim_fft4_avx_fma); }
247#endif
248
249TEST(fft, reim_fft2_ref_vs_naive) { test_reim_fft_ref_vs_naive<2>(fill_reim_fft2_omegas, reim_fft2_ref); }
250
251TEST(fft, reim_fft_bfs_16_ref_vs_naive) {
252 for (const uint64_t m : {16, 32, 64, 128, 256, 512, 1024, 2048}) {
253 std::vector<double> om(2 * m);
254 std::vector<double> data(2 * m);
255 std::vector<double> datacopy(2 * m);
256 double* omg = om.data();
257 fill_reim_fft_bfs_16_omegas(m, 0.25, &omg);
258 ASSERT_LE(omg - om.data(), ptrdiff_t(2 * m)); // it may depend on m
259 for (uint64_t i = 0; i < m; ++i) {
260 datacopy[i] = data[i] = (rand() % 100) - 50;
261 datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
262 }
263 omg = om.data();
264 reim_fft_bfs_16_ref(m, datacopy.data(), datacopy.data() + m, &omg);
265 reim_naive_fft(m, 0.25, data.data(), data.data() + m);
266 double d = 0;
267 for (uint64_t i = 0; i < 2 * m; ++i) {
268 d += fabs(datacopy[i] - data[i]);
269 }
270 ASSERT_LE(d, 1e-7);
271 }
272}
273
274TEST(fft, reim_fft_rec_16_ref_vs_naive) {
275 for (const uint64_t m : {2048, 4096, 8192, 32768, 65536}) {
276 std::vector<double> om(2 * m);
277 std::vector<double> data(2 * m);
278 std::vector<double> datacopy(2 * m);
279 double* omg = om.data();
280 fill_reim_fft_rec_16_omegas(m, 0.25, &omg);
281 ASSERT_LE(omg - om.data(), ptrdiff_t(2 * m)); // it may depend on m
282 for (uint64_t i = 0; i < m; ++i) {
283 datacopy[i] = data[i] = (rand() % 100) - 50;
284 datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
285 }
286 omg = om.data();
287 reim_fft_rec_16_ref(m, datacopy.data(), datacopy.data() + m, &omg);
288 reim_naive_fft(m, 0.25, data.data(), data.data() + m);
289 double d = 0;
290 for (uint64_t i = 0; i < 2 * m; ++i) {
291 d += fabs(datacopy[i] - data[i]);
292 }
293 ASSERT_LE(d, 1e-5);
294 }
295}
296
297TEST(fft, reim_fft_ref_vs_naive) {
298 for (const uint64_t m : {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 32768, 65536}) {
299 std::vector<double> om(2 * m);
300 std::vector<double> data(2 * m);
301 std::vector<double> datacopy(2 * m);
302 REIM_FFT_PRECOMP* precomp = new_reim_fft_precomp(m, 0);
303 for (uint64_t i = 0; i < m; ++i) {
304 datacopy[i] = data[i] = (rand() % 100) - 50;
305 datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
306 }
307 reim_fft_ref(precomp, datacopy.data());
308 reim_naive_fft(m, 0.25, data.data(), data.data() + m);
309 double d = 0;
310 for (uint64_t i = 0; i < 2 * m; ++i) {
311 d += fabs(datacopy[i] - data[i]);
312 }
313 ASSERT_LE(d, 1e-5) << m;
314 delete_reim_fft_precomp(precomp);
315 }
316}
317
318#ifdef __aarch64__
319EXPORT REIM_FFT_PRECOMP* new_reim_fft_precomp_neon(uint32_t m, uint32_t num_buffers);
320EXPORT void reim_fft_neon(const REIM_FFT_PRECOMP* precomp, double* d);
321TEST(fft, reim_fft_neon_vs_naive) {
322 for (const uint64_t m : {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 32768, 65536}) {
323 std::vector<double> om(2 * m);
324 std::vector<double> data(2 * m);
325 std::vector<double> datacopy(2 * m);
326 REIM_FFT_PRECOMP* precomp = new_reim_fft_precomp_neon(m, 0);
327 for (uint64_t i = 0; i < m; ++i) {
328 datacopy[i] = data[i] = (rand() % 100) - 50;
329 datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
330 }
331 reim_fft_neon(precomp, datacopy.data());
332 reim_naive_fft(m, 0.25, data.data(), data.data() + m);
333 double d = 0;
334 for (uint64_t i = 0; i < 2 * m; ++i) {
335 d += fabs(datacopy[i] - data[i]);
336 }
337 ASSERT_LE(d, 1e-5) << m;
338 delete_reim_fft_precomp(precomp);
339 }
340}
341#endif
342
343typedef void (*FILL_REIM_IFFT_OMG_F)(const double entry_pwr, double** omg);
344typedef void (*REIM_IFFT_F)(double* dre, double* dim, const void* omega);
345
346// template to test a fixed-dimension fft vs. naive
347template <uint64_t N>
348void test_reim_ifft_ref_vs_naive(FILL_REIM_IFFT_OMG_F fill_omega_f, REIM_IFFT_F reim_ifft_f) {
349 double om[N];
350 double data[2 * N];
351 double datacopy[2 * N];
352 double* omg = om;
353 fill_omega_f(0.25, &omg);
354 ASSERT_EQ(omg - om, ptrdiff_t(N)); // it may depend on N
355 for (uint64_t i = 0; i < N; ++i) {
356 datacopy[i] = data[i] = (rand() % 100) - 50;
357 datacopy[N + i] = data[N + i] = (rand() % 100) - 50;
358 }
359 reim_ifft_f(datacopy, datacopy + N, om);
360 reim_naive_ifft(N, 0.25, data, data + N);
361 double d = 0;
362 for (uint64_t i = 0; i < 2 * N; ++i) {
363 d += fabs(datacopy[i] - data[i]);
364 }
365 ASSERT_LE(d, 1e-7);
366}
367
368template <uint64_t N>
369void test_reim_ifft_ref_vs_accel(REIM_IFFT_F reim_ifft_ref_f, REIM_IFFT_F reim_ifft_accel_f) {
370 double om[N];
371 double data[2 * N];
372 double datacopy[2 * N];
373 for (uint64_t i = 0; i < N; ++i) {
374 om[i] = (rand() % 100) - 50;
375 datacopy[i] = data[i] = (rand() % 100) - 50;
376 datacopy[N + i] = data[N + i] = (rand() % 100) - 50;
377 }
378 reim_ifft_ref_f(datacopy, datacopy + N, om);
379 reim_ifft_accel_f(data, data + N, om);
380 double d = 0;
381 for (uint64_t i = 0; i < 2 * N; ++i) {
382 d += fabs(datacopy[i] - data[i]);
383 }
384 if (d > 1e-15) {
385 for (uint64_t i = 0; i < N; ++i) {
386 printf("%" PRId64 " %lf %lf %lf %lf\n", i, data[i], data[N + i], datacopy[i], datacopy[N + i]);
387 }
388 ASSERT_LE(d, 0);
389 }
390}
391
392TEST(fft, reim_ifft16_ref_vs_naive) { test_reim_ifft_ref_vs_naive<16>(fill_reim_ifft16_omegas, reim_ifft16_ref); }
393
394#ifdef __x86_64__
395TEST(fft, reim_ifft16_ref_vs_fma) { test_reim_ifft_ref_vs_accel<16>(reim_ifft16_ref, reim_ifft16_avx_fma); }
396#endif
397
398TEST(fft, reim_ifft8_ref_vs_naive) { test_reim_ifft_ref_vs_naive<8>(fill_reim_ifft8_omegas, reim_ifft8_ref); }
399
400#ifdef __x86_64__
401TEST(fft, reim_ifft8_ref_vs_fma) { test_reim_ifft_ref_vs_accel<8>(reim_ifft8_ref, reim_ifft8_avx_fma); }
402#endif
403
404TEST(fft, reim_ifft4_ref_vs_naive) { test_reim_ifft_ref_vs_naive<4>(fill_reim_ifft4_omegas, reim_ifft4_ref); }
405
406#ifdef __x86_64__
407TEST(fft, reim_ifft4_ref_vs_fma) { test_reim_ifft_ref_vs_accel<4>(reim_ifft4_ref, reim_ifft4_avx_fma); }
408#endif
409
410TEST(fft, reim_ifft2_ref_vs_naive) { test_reim_ifft_ref_vs_naive<2>(fill_reim_ifft2_omegas, reim_ifft2_ref); }
411
412TEST(fft, reim_ifft_bfs_16_ref_vs_naive) {
413 for (const uint64_t m : {16, 32, 64, 128, 256, 512, 1024, 2048}) {
414 std::vector<double> om(2 * m);
415 std::vector<double> data(2 * m);
416 std::vector<double> datacopy(2 * m);
417 double* omg = om.data();
418 fill_reim_ifft_bfs_16_omegas(m, 0.25, &omg);
419 ASSERT_LE(omg - om.data(), ptrdiff_t(2 * m)); // it may depend on m
420 for (uint64_t i = 0; i < m; ++i) {
421 datacopy[i] = data[i] = (rand() % 100) - 50;
422 datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
423 }
424 omg = om.data();
425 reim_ifft_bfs_16_ref(m, datacopy.data(), datacopy.data() + m, &omg);
426 reim_naive_ifft(m, 0.25, data.data(), data.data() + m);
427 double d = 0;
428 for (uint64_t i = 0; i < 2 * m; ++i) {
429 d += fabs(datacopy[i] - data[i]);
430 }
431 ASSERT_LE(d, 1e-7);
432 }
433}
434
435TEST(fft, reim_ifft_rec_16_ref_vs_naive) {
436 for (const uint64_t m : {2048, 4096, 8192, 32768, 65536}) {
437 std::vector<double> om(2 * m);
438 std::vector<double> data(2 * m);
439 std::vector<double> datacopy(2 * m);
440 double* omg = om.data();
441 fill_reim_ifft_rec_16_omegas(m, 0.25, &omg);
442 ASSERT_LE(omg - om.data(), ptrdiff_t(2 * m)); // it may depend on m
443 for (uint64_t i = 0; i < m; ++i) {
444 datacopy[i] = data[i] = (rand() % 100) - 50;
445 datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
446 }
447 omg = om.data();
448 reim_ifft_rec_16_ref(m, datacopy.data(), datacopy.data() + m, &omg);
449 reim_naive_ifft(m, 0.25, data.data(), data.data() + m);
450 double d = 0;
451 for (uint64_t i = 0; i < 2 * m; ++i) {
452 d += fabs(datacopy[i] - data[i]);
453 }
454 ASSERT_LE(d, 1e-5);
455 }
456}
457
458TEST(fft, reim_ifft_ref_vs_naive) {
459 for (const uint64_t m : {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 32768, 65536}) {
460 std::vector<double> om(2 * m);
461 std::vector<double> data(2 * m);
462 std::vector<double> datacopy(2 * m);
463 REIM_IFFT_PRECOMP* precomp = new_reim_ifft_precomp(m, 0);
464 for (uint64_t i = 0; i < m; ++i) {
465 datacopy[i] = data[i] = (rand() % 100) - 50;
466 datacopy[m + i] = data[m + i] = (rand() % 100) - 50;
467 }
468 reim_ifft_ref(precomp, datacopy.data());
469 reim_naive_ifft(m, 0.25, data.data(), data.data() + m);
470 double d = 0;
471 for (uint64_t i = 0; i < 2 * m; ++i) {
472 d += fabs(datacopy[i] - data[i]);
473 }
474 ASSERT_LE(d, 1e-5) << m;
475 delete_reim_ifft_precomp(precomp);
476 }
477}