poulpy/spqlios/lib/test/spqlios_cplx_test.cpp
spqlios/lib/test/spqlios_cplx_test.cppRaw
1#include <cmath>
2
3#include "gtest/gtest.h"
4#include "spqlios/commons_private.h"
5#include "spqlios/cplx/cplx_fft.h"
6#include "spqlios/cplx/cplx_fft_internal.h"
7#include "spqlios/cplx/cplx_fft_private.h"
8
9#ifdef __x86_64__
10TEST(fft, ifft16_fma_vs_ref) {
11 CPLX data[16];
12 CPLX omega[8];
13 for (uint64_t i = 0; i < 32; ++i) ((double*)data)[i] = 2 * i + 1; //(rand()%100)-50;
14 for (uint64_t i = 0; i < 16; ++i) ((double*)omega)[i] = i + 1; //(rand()%100)-50;
15 CPLX copydata[16];
16 CPLX copyomega[8];
17 memcpy(copydata, data, sizeof(copydata));
18 memcpy(copyomega, omega, sizeof(copyomega));
19 cplx_ifft16_avx_fma(data, omega);
20 cplx_ifft16_ref(copydata, copyomega);
21 double distance = 0;
22 for (uint64_t i = 0; i < 16; ++i) {
23 double d1 = fabs(data[i][0] - copydata[i][0]);
24 double d2 = fabs(data[i][0] - copydata[i][0]);
25 if (d1 > distance) distance = d1;
26 if (d2 > distance) distance = d2;
27 }
28 /*
29 printf("data:\n");
30 for (uint64_t i=0; i<4; ++i) {
31 for (uint64_t j=0; j<8; ++j) {
32 printf("%.5lf ", data[4 * i + j / 2][j % 2]);
33 }
34 printf("\n");
35 }
36 printf("copydata:\n");
37 for (uint64_t i=0; i<4; ++i) {
38 for (uint64_t j=0; j<8; ++j) {
39 printf("%5.5lf ", copydata[4 * i + j / 2][j % 2]);
40 }
41 printf("\n");
42 }
43 */
44 ASSERT_EQ(distance, 0);
45}
46
47#endif
48
49void cplx_zero(CPLX r) { r[0] = r[1] = 0; }
50void cplx_addmul(CPLX r, const CPLX a, const CPLX b) {
51 double re = r[0] + a[0] * b[0] - a[1] * b[1];
52 double im = r[1] + a[0] * b[1] + a[1] * b[0];
53 r[0] = re;
54 r[1] = im;
55}
56
57void halfcfft_eval(CPLX res, uint32_t nn, uint32_t k, const CPLX* coeffs, const CPLX* powomegas) {
58 const uint32_t N = nn / 2;
59 cplx_zero(res);
60 for (uint64_t i = 0; i < N; ++i) {
61 cplx_addmul(res, coeffs[i], powomegas[(k * i) % (2 * nn)]);
62 }
63}
64void halfcfft_naive(uint32_t nn, CPLX* data) {
65 const uint32_t N = nn / 2;
66 CPLX* in = (CPLX*)malloc(N * sizeof(CPLX));
67 CPLX* powomega = (CPLX*)malloc(2 * nn * sizeof(CPLX));
68 for (uint64_t i = 0; i < (2 * nn); ++i) {
69 powomega[i][0] = m_accurate_cos((M_PI * i) / nn);
70 powomega[i][1] = m_accurate_sin((M_PI * i) / nn);
71 }
72 memcpy(in, data, N * sizeof(CPLX));
73 for (uint64_t j = 0; j < N; ++j) {
74 uint64_t p = rint(log2(N)) + 2;
75 uint64_t k = revbits(p, j) + 1;
76 halfcfft_eval(data[j], nn, k, in, powomega);
77 }
78 free(powomega);
79 free(in);
80}
81
82#ifdef __x86_64__
83TEST(fft, fft16_fma_vs_ref) {
84 CPLX data[16];
85 CPLX omega[8];
86 for (uint64_t i = 0; i < 32; ++i) ((double*)data)[i] = rand() % 1000;
87 for (uint64_t i = 0; i < 16; ++i) ((double*)omega)[i] = rand() % 1000;
88 CPLX copydata[16];
89 CPLX copyomega[8];
90 memcpy(copydata, data, sizeof(copydata));
91 memcpy(copyomega, omega, sizeof(copyomega));
92 cplx_fft16_avx_fma(data, omega);
93 cplx_fft16_ref(copydata, omega);
94 double distance = 0;
95 for (uint64_t i = 0; i < 16; ++i) {
96 double d1 = fabs(data[i][0] - copydata[i][0]);
97 double d2 = fabs(data[i][0] - copydata[i][0]);
98 if (d1 > distance) distance = d1;
99 if (d2 > distance) distance = d2;
100 }
101 ASSERT_EQ(distance, 0);
102}
103#endif
104
105TEST(fft, citwiddle_then_invcitwiddle) {
106 CPLX om;
107 CPLX ombar;
108 CPLX data[2];
109 CPLX copydata[2];
110 om[0] = cos(3);
111 om[1] = sin(3);
112 ombar[0] = om[0];
113 ombar[1] = -om[1];
114 data[0][0] = 47;
115 data[0][1] = 23;
116 data[1][0] = -12;
117 data[1][1] = -9;
118 memcpy(copydata, data, sizeof(copydata));
119 citwiddle(data[0], data[1], om);
120 invcitwiddle(data[0], data[1], ombar);
121 double distance = 0;
122 for (uint64_t i = 0; i < 2; ++i) {
123 double d1 = fabs(data[i][0] - 2 * copydata[i][0]);
124 double d2 = fabs(data[i][1] - 2 * copydata[i][1]);
125 if (d1 > distance) distance = d1;
126 if (d2 > distance) distance = d2;
127 }
128 ASSERT_LE(distance, 1e-9);
129}
130
131TEST(fft, ctwiddle_then_invctwiddle) {
132 CPLX om;
133 CPLX ombar;
134 CPLX data[2];
135 CPLX copydata[2];
136 om[0] = cos(3);
137 om[1] = sin(3);
138 ombar[0] = om[0];
139 ombar[1] = -om[1];
140 data[0][0] = 47;
141 data[0][1] = 23;
142 data[1][0] = -12;
143 data[1][1] = -9;
144 memcpy(copydata, data, sizeof(copydata));
145 ctwiddle(data[0], data[1], om);
146 invctwiddle(data[0], data[1], ombar);
147 double distance = 0;
148 for (uint64_t i = 0; i < 2; ++i) {
149 double d1 = fabs(data[i][0] - 2 * copydata[i][0]);
150 double d2 = fabs(data[i][1] - 2 * copydata[i][1]);
151 if (d1 > distance) distance = d1;
152 if (d2 > distance) distance = d2;
153 }
154 ASSERT_LE(distance, 1e-9);
155}
156
157TEST(fft, fft16_then_ifft16_ref) {
158 CPLX full_omegas[64];
159 CPLX full_omegabars[64];
160 for (uint64_t i = 0; i < 64; ++i) {
161 full_omegas[i][0] = cos(M_PI * i / 32.);
162 full_omegas[i][1] = sin(M_PI * i / 32.);
163 full_omegabars[i][0] = full_omegas[i][0];
164 full_omegabars[i][1] = -full_omegas[i][1];
165 }
166 CPLX omega[8];
167 CPLX omegabar[8];
168 cplx_set(omega[0], full_omegas[8]); // j
169 cplx_set(omega[1], full_omegas[4]); // k
170 cplx_set(omega[2], full_omegas[2]); // l
171 cplx_set(omega[3], full_omegas[10]); // lj
172 cplx_set(omega[4], full_omegas[1]); // n
173 cplx_set(omega[5], full_omegas[9]); // nj
174 cplx_set(omega[6], full_omegas[5]); // nk
175 cplx_set(omega[7], full_omegas[13]); // njk
176 cplx_set(omegabar[0], full_omegabars[1]); // n
177 cplx_set(omegabar[1], full_omegabars[9]); // nj
178 cplx_set(omegabar[2], full_omegabars[5]); // nk
179 cplx_set(omegabar[3], full_omegabars[13]); // njk
180 cplx_set(omegabar[4], full_omegabars[2]); // l
181 cplx_set(omegabar[5], full_omegabars[10]); // lj
182 cplx_set(omegabar[6], full_omegabars[4]); // k
183 cplx_set(omegabar[7], full_omegabars[8]); // j
184 CPLX data[16];
185 CPLX copydata[16];
186 for (uint64_t i = 0; i < 32; ++i) ((double*)data)[i] = rand() % 1000;
187 memcpy(copydata, data, sizeof(copydata));
188 cplx_fft16_ref(data, omega);
189 cplx_ifft16_ref(data, omegabar);
190 double distance = 0;
191 for (uint64_t i = 0; i < 16; ++i) {
192 double d1 = fabs(data[i][0] - 16 * copydata[i][0]);
193 double d2 = fabs(data[i][0] - 16 * copydata[i][0]);
194 if (d1 > distance) distance = d1;
195 if (d2 > distance) distance = d2;
196 }
197 ASSERT_LE(distance, 1e-9);
198}
199
200TEST(fft, halfcfft_ref_vs_naive) {
201 for (uint64_t nn : {4, 8, 16, 64, 256, 8192}) {
202 uint64_t m = nn / 2;
203 CPLX_FFT_PRECOMP* tables = new_cplx_fft_precomp(m, 0);
204 CPLX* a = (CPLX*)spqlios_alloc_custom_align(32, m * sizeof(CPLX));
205 CPLX* a1 = (CPLX*)spqlios_alloc_custom_align(32, m * sizeof(CPLX));
206 CPLX* a2 = (CPLX*)spqlios_alloc_custom_align(32, m * sizeof(CPLX));
207 int64_t p = 1 << 16;
208 for (uint32_t i = 0; i < m; i++) {
209 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
210 a[i][1] = (rand() % p) - p / 2;
211 }
212 memcpy(a1, a, m * sizeof(CPLX));
213 memcpy(a2, a, m * sizeof(CPLX));
214
215 halfcfft_naive(nn, a1);
216 cplx_fft_naive(m, 0.25, a2);
217
218 double d = 0;
219 for (uint32_t i = 0; i < nn / 2; i++) {
220 double dre = fabs(a1[i][0] - a2[i][0]);
221 double dim = fabs(a1[i][1] - a2[i][1]);
222 if (dre > d) d = dre;
223 if (dim > d) d = dim;
224 }
225 ASSERT_LE(d, nn * 1e-10) << nn;
226 spqlios_free(a);
227 spqlios_free(a1);
228 spqlios_free(a2);
229 delete_cplx_fft_precomp(tables);
230 }
231}
232
233#ifdef __x86_64__
234TEST(fft, halfcfft_fma_vs_ref) {
235 typedef void (*FFTF)(const CPLX_FFT_PRECOMP*, void* data);
236 for (FFTF fft : {cplx_fft_ref, cplx_fft_avx2_fma}) {
237 for (uint64_t nn : {8, 16, 32, 64, 1024, 8192, 65536}) {
238 uint64_t m = nn / 2;
239 CPLX_FFT_PRECOMP* tables = new_cplx_fft_precomp(m, 0);
240 CPLX* a = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
241 CPLX* a1 = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
242 CPLX* a2 = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
243 int64_t p = 1 << 16;
244 for (uint32_t i = 0; i < nn / 2; i++) {
245 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
246 a[i][1] = (rand() % p) - p / 2;
247 }
248 memcpy(a1, a, nn / 2 * sizeof(CPLX));
249 memcpy(a2, a, nn / 2 * sizeof(CPLX));
250 cplx_fft_naive(m, 0.25, a2);
251 fft(tables, a1);
252 double d = 0;
253 for (uint32_t i = 0; i < nn / 2; i++) {
254 double dre = fabs(a1[i][0] - a2[i][0]);
255 double dim = fabs(a1[i][1] - a2[i][1]);
256 if (dre > d) d = dre;
257 if (dim > d) d = dim;
258 }
259 ASSERT_LE(d, nn * 1e-10) << nn;
260 spqlios_free(a);
261 spqlios_free(a1);
262 spqlios_free(a2);
263 delete_cplx_fft_precomp(tables);
264 }
265 }
266}
267#endif
268
269TEST(fft, halfcfft_then_ifft_ref) {
270 for (uint64_t nn : {4, 8, 16, 32, 64, 1024, 8192, 65536}) {
271 uint64_t m = nn / 2;
272 CPLX_FFT_PRECOMP* tables = new_cplx_fft_precomp(m, 0);
273 CPLX_IFFT_PRECOMP* itables = new_cplx_ifft_precomp(m, 0);
274 CPLX* a = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
275 CPLX* a1 = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
276 int64_t p = 1 << 16;
277 for (uint32_t i = 0; i < nn / 2; i++) {
278 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
279 a[i][1] = (rand() % p) - p / 2;
280 }
281 memcpy(a1, a, nn / 2 * sizeof(CPLX));
282 cplx_fft_ref(tables, a1);
283 cplx_ifft_ref(itables, a1);
284 double d = 0;
285 for (uint32_t i = 0; i < nn / 2; i++) {
286 double dre = fabs(a[i][0] - a1[i][0] / (nn / 2));
287 double dim = fabs(a[i][1] - a1[i][1] / (nn / 2));
288 if (dre > d) d = dre;
289 if (dim > d) d = dim;
290 }
291 ASSERT_LE(d, 1e-8);
292 spqlios_free(a);
293 spqlios_free(a1);
294 delete_cplx_fft_precomp(tables);
295 delete_cplx_ifft_precomp(itables);
296 }
297}
298
299#ifdef __x86_64__
300TEST(fft, halfcfft_ifft_fma_vs_ref) {
301 for (IFFT_FUNCTION ifft : {cplx_ifft_ref, cplx_ifft_avx2_fma}) {
302 for (uint64_t nn : {8, 16, 32, 1024, 4096, 8192, 65536}) {
303 uint64_t m = nn / 2;
304 CPLX_IFFT_PRECOMP* itables = new_cplx_ifft_precomp(m, 0);
305 CPLX* a = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
306 CPLX* a1 = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
307 CPLX* a2 = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
308 int64_t p = 1 << 16;
309 for (uint32_t i = 0; i < nn / 2; i++) {
310 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
311 a[i][1] = (rand() % p) - p / 2;
312 }
313 memcpy(a1, a, nn / 2 * sizeof(CPLX));
314 memcpy(a2, a, nn / 2 * sizeof(CPLX));
315 cplx_ifft_naive(m, 0.25, a2);
316 ifft(itables, a1);
317 double d = 0;
318 for (uint32_t i = 0; i < nn / 2; i++) {
319 double dre = fabs(a1[i][0] - a2[i][0]);
320 double dim = fabs(a1[i][1] - a2[i][1]);
321 if (dre > d) d = dre;
322 if (dim > d) d = dim;
323 }
324 ASSERT_LE(d, 1e-8);
325 spqlios_free(a);
326 spqlios_free(a1);
327 spqlios_free(a2);
328 delete_cplx_ifft_precomp(itables);
329 }
330 }
331}
332#endif
333
334// test the reference and simple implementations of mul on all dimensions
335TEST(fftvec, cplx_fftvec_mul_ref) {
336 for (uint64_t nn : {2, 4, 8, 16, 32, 1024, 4096, 8192, 65536}) {
337 uint64_t m = nn / 2;
338 CPLX_FFTVEC_MUL_PRECOMP* precomp = new_cplx_fftvec_mul_precomp(m);
339 CPLX* a = new CPLX[m];
340 CPLX* b = new CPLX[m];
341 CPLX* r0 = new CPLX[m];
342 CPLX* r1 = new CPLX[m];
343 CPLX* r2 = new CPLX[m];
344 int64_t p = 1 << 16;
345 for (uint32_t i = 0; i < m; i++) {
346 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
347 a[i][1] = (rand() % p) - p / 2;
348 b[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
349 b[i][1] = (rand() % p) - p / 2;
350 r2[i][0] = r1[i][0] = r0[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
351 r2[i][1] = r1[i][1] = r0[i][1] = (rand() % p) - p / 2;
352 }
353 cplx_fftvec_mul_simple(m, r0, a, b);
354 cplx_fftvec_mul_ref(precomp, r1, a, b);
355 for (uint32_t i = 0; i < m; i++) {
356 r2[i][0] = a[i][0] * b[i][0] - a[i][1] * b[i][1];
357 r2[i][1] = a[i][0] * b[i][1] + a[i][1] * b[i][0];
358 ASSERT_LE(fabs(r1[i][0] - r2[i][0]) + fabs(r1[i][1] - r2[i][1]), 1e-8);
359 ASSERT_LE(fabs(r0[i][0] - r2[i][0]) + fabs(r0[i][1] - r2[i][1]), 1e-8);
360 }
361 delete[] a;
362 delete[] b;
363 delete[] r0;
364 delete[] r1;
365 delete[] r2;
366 delete_cplx_fftvec_mul_precomp(precomp);
367 }
368}
369
370// test the reference and simple implementations of addmul on all dimensions
371TEST(fftvec, cplx_fftvec_addmul_ref) {
372 for (uint64_t nn : {2, 4, 8, 16, 32, 1024, 4096, 8192, 65536}) {
373 uint64_t m = nn / 2;
374 CPLX_FFTVEC_ADDMUL_PRECOMP* precomp = new_cplx_fftvec_addmul_precomp(m);
375 CPLX* a = new CPLX[m];
376 CPLX* b = new CPLX[m];
377 CPLX* r0 = new CPLX[m];
378 CPLX* r1 = new CPLX[m];
379 CPLX* r2 = new CPLX[m];
380 int64_t p = 1 << 16;
381 for (uint32_t i = 0; i < m; i++) {
382 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
383 a[i][1] = (rand() % p) - p / 2;
384 b[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
385 b[i][1] = (rand() % p) - p / 2;
386 r2[i][0] = r1[i][0] = r0[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
387 r2[i][1] = r1[i][1] = r0[i][1] = (rand() % p) - p / 2;
388 }
389 cplx_fftvec_addmul_simple(m, r0, a, b);
390 cplx_fftvec_addmul_ref(precomp, r1, a, b);
391 for (uint32_t i = 0; i < m; i++) {
392 r2[i][0] += a[i][0] * b[i][0] - a[i][1] * b[i][1];
393 r2[i][1] += a[i][0] * b[i][1] + a[i][1] * b[i][0];
394 ASSERT_LE(fabs(r1[i][0] - r2[i][0]) + fabs(r1[i][1] - r2[i][1]), 1e-8);
395 ASSERT_LE(fabs(r0[i][0] - r2[i][0]) + fabs(r0[i][1] - r2[i][1]), 1e-8);
396 }
397 delete[] a;
398 delete[] b;
399 delete[] r0;
400 delete[] r1;
401 delete[] r2;
402 delete_cplx_fftvec_addmul_precomp(precomp);
403 }
404}
405
406// comparative tests between mul ref vs. optimized (only relevant dimensions)
407TEST(fftvec, cplx_fftvec_mul_ref_vs_optim) {
408 struct totest {
409 FFTVEC_MUL_FUNCTION f;
410 uint64_t min_m;
411 totest(FFTVEC_MUL_FUNCTION f, uint64_t min_m) : f(f), min_m(min_m) {}
412 };
413 std::vector<totest> totestset;
414 totestset.emplace_back(cplx_fftvec_mul, 1);
415#ifdef __x86_64__
416 totestset.emplace_back(cplx_fftvec_mul_fma, 8);
417#endif
418 for (uint64_t m : {1, 2, 4, 8, 16, 1024, 4096, 8192, 65536}) {
419 CPLX_FFTVEC_MUL_PRECOMP* precomp = new_cplx_fftvec_mul_precomp(m);
420 for (const totest& t : totestset) {
421 if (t.min_m > m) continue;
422 CPLX* a = new CPLX[m];
423 CPLX* b = new CPLX[m];
424 CPLX* r1 = new CPLX[m];
425 CPLX* r2 = new CPLX[m];
426 int64_t p = 1 << 16;
427 for (uint32_t i = 0; i < m; i++) {
428 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
429 a[i][1] = (rand() % p) - p / 2;
430 b[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
431 b[i][1] = (rand() % p) - p / 2;
432 r2[i][0] = r1[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
433 r2[i][1] = r1[i][1] = (rand() % p) - p / 2;
434 }
435 t.f(precomp, r1, a, b);
436 cplx_fftvec_mul_ref(precomp, r2, a, b);
437 for (uint32_t i = 0; i < m; i++) {
438 double dre = fabs(r1[i][0] - r2[i][0]);
439 double dim = fabs(r1[i][1] - r2[i][1]);
440 ASSERT_LE(dre, 1e-8);
441 ASSERT_LE(dim, 1e-8);
442 }
443 delete[] a;
444 delete[] b;
445 delete[] r1;
446 delete[] r2;
447 }
448 delete_cplx_fftvec_mul_precomp(precomp);
449 }
450}
451
452// comparative tests between addmul ref vs. optimized (only relevant dimensions)
453TEST(fftvec, cplx_fftvec_addmul_ref_vs_optim) {
454 struct totest {
455 FFTVEC_ADDMUL_FUNCTION f;
456 uint64_t min_m;
457 totest(FFTVEC_ADDMUL_FUNCTION f, uint64_t min_m) : f(f), min_m(min_m) {}
458 };
459 std::vector<totest> totestset;
460 totestset.emplace_back(cplx_fftvec_addmul, 1);
461#ifdef __x86_64__
462 totestset.emplace_back(cplx_fftvec_addmul_fma, 8);
463#endif
464 for (uint64_t m : {1, 2, 4, 8, 16, 1024, 4096, 8192, 65536}) {
465 CPLX_FFTVEC_ADDMUL_PRECOMP* precomp = new_cplx_fftvec_addmul_precomp(m);
466 for (const totest& t : totestset) {
467 if (t.min_m > m) continue;
468 CPLX* a = new CPLX[m];
469 CPLX* b = new CPLX[m];
470 CPLX* r1 = new CPLX[m];
471 CPLX* r2 = new CPLX[m];
472 int64_t p = 1 << 16;
473 for (uint32_t i = 0; i < m; i++) {
474 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
475 a[i][1] = (rand() % p) - p / 2;
476 b[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
477 b[i][1] = (rand() % p) - p / 2;
478 r2[i][0] = r1[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
479 r2[i][1] = r1[i][1] = (rand() % p) - p / 2;
480 }
481 t.f(precomp, r1, a, b);
482 cplx_fftvec_addmul_ref(precomp, r2, a, b);
483 for (uint32_t i = 0; i < m; i++) {
484 double dre = fabs(r1[i][0] - r2[i][0]);
485 double dim = fabs(r1[i][1] - r2[i][1]);
486 ASSERT_LE(dre, 1e-8);
487 ASSERT_LE(dim, 1e-8);
488 }
489 delete[] a;
490 delete[] b;
491 delete[] r1;
492 delete[] r2;
493 }
494 delete_cplx_fftvec_addmul_precomp(precomp);
495 }
496}