arnaucube
arnaucube/poulpyPublic Self-hosted repository
poulpy/poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/test/spqlios_cplx_test.cpp
poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/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}