poulpy/spqlios/lib/test/spqlios_test.cpp
spqlios/lib/test/spqlios_test.cppRaw
1#include <cassert>
2#include <cinttypes>
3#include <cmath>
4#include <complex>
5#include <cstdint>
6#include <cstdlib>
7#include <iostream>
8#include <string>
9#include <vector>
10
11#include "gtest/gtest.h"
12#include "spqlios/cplx/cplx_fft_internal.h"
13
14using namespace std;
15
16/*namespace {
17bool very_close(const double& a, const double& b) {
18 bool reps = (abs(a - b) < 1e-5);
19 if (!reps) {
20 cerr << "not close: " << a << " vs. " << b << endl;
21 }
22 return reps;
23}
24
25}*/ // namespace
26
27TEST(fft, fftvec_convolution) {
28 uint64_t nn = 65536; // vary accross (8192, 16384), 32768, 65536
29 static const uint64_t k = 18; // vary from 10 to 20
30 // double* buf_fft = fft_precomp_get_buffer(tables, 0);
31 // double* buf_ifft = ifft_precomp_get_buffer(itables, 0);
32 double* a = (double*)spqlios_alloc_custom_align(32, nn * 8);
33 double* a2 = (double*)spqlios_alloc_custom_align(32, nn * 8);
34 double* b = (double*)spqlios_alloc_custom_align(32, nn * 8);
35 double* dist_vector = (double*)spqlios_alloc_custom_align(32, nn * 8);
36 int64_t p = UINT64_C(1) << k;
37 printf("p size: %" PRId64 "\n", p);
38 for (uint32_t i = 0; i < nn; i++) {
39 a[i] = (rand() % p) - p / 2; // between -p/2 and p/2
40 b[i] = (rand() % p) - p / 2;
41 a2[i] = 0;
42 }
43 cplx_fft_simple(nn / 2, a);
44 cplx_fft_simple(nn / 2, b);
45 cplx_fftvec_addmul_simple(nn / 2, a2, a, b);
46 cplx_ifft_simple(nn / 2, a2); // normalization is missing
47 double distance = 0;
48 // for (int32_t i = 0; i < 10; i++) {
49 // printf("%lf %lf\n", a2[i], a2[i] / (nn / 2.));
50 //}
51 for (uint32_t i = 0; i < nn; i++) {
52 double curdist = fabs(a2[i] / (nn / 2.) - rint(a2[i] / (nn / 2.)));
53 if (distance < curdist) distance = curdist;
54 dist_vector[i] = a2[i] / (nn / 2.) - rint(a2[i] / (nn / 2.));
55 }
56 printf("distance: %lf\n", distance);
57 ASSERT_LE(distance, 0.5); // switch from previous 0.1 to 0.5 per experiment 1 reqs
58 // double a3[] = {2,4,4,4,5,5,7,9}; //instead of dist_vector, for test
59 // nn = 8;
60 double mean = 0;
61 for (uint32_t i = 0; i < nn; i++) {
62 mean = mean + dist_vector[i];
63 }
64 mean = mean / nn;
65 double variance = 0;
66 for (uint32_t i = 0; i < nn; i++) {
67 variance = variance + pow((mean - dist_vector[i]), 2);
68 }
69 double stdev = sqrt(variance / nn);
70 printf("stdev: %lf\n", stdev);
71
72 spqlios_free(a);
73 spqlios_free(b);
74 spqlios_free(a2);
75 spqlios_free(dist_vector);
76}
77
78typedef double CPLX[2];
79EXPORT uint32_t revbits(uint32_t i, uint32_t v);
80
81void cplx_zero(CPLX r);
82void cplx_addmul(CPLX r, const CPLX a, const CPLX b);
83
84void halfcfft_eval(CPLX res, uint32_t nn, uint32_t k, const CPLX* coeffs, const CPLX* powomegas);
85void halfcfft_naive(uint32_t nn, CPLX* data);
86
87EXPORT void cplx_set(CPLX, const CPLX);
88EXPORT void citwiddle(CPLX, CPLX, const CPLX);
89EXPORT void invcitwiddle(CPLX, CPLX, const CPLX);
90EXPORT void ctwiddle(CPLX, CPLX, const CPLX);
91EXPORT void invctwiddle(CPLX, CPLX, const CPLX);
92
93#include "../spqlios/cplx/cplx_fft_private.h"
94#include "../spqlios/reim/reim_fft_internal.h"
95#include "../spqlios/reim/reim_fft_private.h"
96#include "../spqlios/reim4/reim4_fftvec_internal.h"
97#include "../spqlios/reim4/reim4_fftvec_private.h"
98
99TEST(fft, simple_fft_test) { // test for checking the simple_fft api
100 uint64_t nn = 8; // vary accross (8192, 16384), 32768, 65536
101 // double* buf_fft = fft_precomp_get_buffer(tables, 0);
102 // double* buf_ifft = ifft_precomp_get_buffer(itables, 0);
103
104 // define the complex coefficients of two polynomials mod X^4-i
105 double a[4][2] = {{1.1, 2.2}, {3.3, 4.4}, {5.5, 6.6}, {7.7, 8.8}};
106 double b[4][2] = {{9., 10.}, {11., 12.}, {13., 14.}, {15., 16.}};
107 double c[4][2]; // for the result
108 double a2[4][2]; // for testing inverse fft
109 memcpy(a2, a, 8 * nn);
110 cplx_fft_simple(4, a);
111 cplx_fft_simple(4, b);
112 cplx_fftvec_mul_simple(4, c, a, b);
113 cplx_ifft_simple(4, c);
114 // c contains the complex coefficients 4.a*b mod X^4-i
115 cplx_ifft_simple(4, a);
116
117 double distance = 0;
118 for (uint32_t i = 0; i < nn / 2; i++) {
119 double dist = fabs(a[i][0] / 4. - a2[i][0]);
120 if (distance < dist) distance = dist;
121 dist = fabs(a[i][1] / 4. - a2[i][1]);
122 if (distance < dist) distance = dist;
123 }
124 printf("distance: %lf\n", distance);
125 ASSERT_LE(distance, 0.1); // switch from previous 0.1 to 0.5 per experiment 1 reqs
126
127 for (uint32_t i = 0; i < nn / 4; i++) {
128 printf("%lf %lf\n", a2[i][0], a[i][0] / (nn / 2.));
129 printf("%lf %lf\n", a2[i][1], a[i][1] / (nn / 2.));
130 }
131}
132
133TEST(fft, reim_test) {
134 // double a[16] __attribute__ ((aligned(32)))= {1.1,2.2,3.3,4.4,5.5,6.6,7.7,8.8,9.9,10.,11.,12.,13.,14.,15.,16.};
135 // double b[16] __attribute__ ((aligned(32)))= {17.,18.,19.,20.,21.,22.,23.,24.,25.,26.,27.,28.,29.,30., 31.,32.};
136 // double c[16] __attribute__ ((aligned(32))); // for the result in reference layout
137 double a[16] = {1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9, 10., 11., 12., 13., 14., 15., 16.};
138 double b[16] = {17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., 30., 31., 32.};
139 double c[16]; // for the result in reference layout
140 reim_fft_simple(8, a);
141 reim_fft_simple(8, b);
142 reim_fftvec_mul_simple(8, c, a, b);
143 reim_ifft_simple(8, c);
144}
145
146TEST(fft, reim_vs_regular_layout_mul_test) {
147 uint64_t nn = 16;
148
149 // define the complex coefficients of two polynomials mod X^8-i
150
151 double a1[8][2] __attribute__((aligned(32))) = {{1.1, 2.2}, {3.3, 4.4}, {5.5, 6.6}, {7.7, 8.8},
152 {9.9, 10.}, {11., 12.}, {13., 14.}, {15., 16.}};
153 double b1[8][2] __attribute__((aligned(32))) = {{17., 18.}, {19., 20.}, {21., 22.}, {23., 24.},
154 {25., 26.}, {27., 28.}, {29., 30.}, {31., 32.}};
155 double c1[8][2] __attribute__((aligned(32))); // for the result
156 double c2[16] __attribute__((aligned(32))); // for the result
157 double c3[8][2] __attribute__((aligned(32))); // for the result
158
159 double* a2 =
160 (double*)spqlios_alloc_custom_align(32, nn / 2 * 2 * sizeof(double)); // for storing the coefs in reim layout
161 double* b2 =
162 (double*)spqlios_alloc_custom_align(32, nn / 2 * 2 * sizeof(double)); // for storing the coefs in reim layout
163 // double* c2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX)); // for storing the coefs in reim
164 // layout
165
166 // organise the coefficients in the reim layout
167 for (uint32_t i = 0; i < nn / 2; i++) {
168 a2[i] = a1[i][0]; // a1 = a2, b1 = b2
169 a2[nn / 2 + i] = a1[i][1];
170 b2[i] = b1[i][0];
171 b2[nn / 2 + i] = b1[i][1];
172 }
173
174 // fft
175 cplx_fft_simple(8, a1);
176 reim_fft_simple(8, a2);
177
178 cplx_fft_simple(8, b1);
179 reim_fft_simple(8, b2);
180
181 cplx_fftvec_mul_simple(8, c1, a1, b1);
182 reim_fftvec_mul_simple(8, c2, a2, b2);
183
184 cplx_ifft_simple(8, c1);
185 reim_ifft_simple(8, c2);
186
187 // check base layout and reim layout result in the same values
188 double d = 0;
189 for (uint32_t i = 0; i < nn / 2; i++) {
190 // printf("RE: cplx_result %lf and reim_result %lf \n", c1[i][0], c2[i]);
191 // printf("IM: cplx_result %lf and reim_result %lf \n", c1[i][1], c2[nn / 2 + i]);
192 double dre = fabs(c1[i][0] - c2[i]);
193 double dim = fabs(c1[i][1] - c2[nn / 2 + i]);
194 if (dre > d) d = dre;
195 if (dim > d) d = dim;
196 ASSERT_LE(d, 1e-7);
197 }
198 ASSERT_LE(d, 1e-7);
199
200 // check converting back to base layout:
201
202 for (uint32_t i = 0; i < nn / 2; i++) {
203 c3[i][0] = c2[i];
204 c3[i][1] = c2[8 + i];
205 }
206
207 d = 0;
208 for (uint32_t i = 0; i < nn / 2; i++) {
209 double dre = fabs(c1[i][0] - c3[i][0]);
210 double dim = fabs(c1[i][1] - c3[i][1]);
211 if (dre > d) d = dre;
212 if (dim > d) d = dim;
213 ASSERT_LE(d, 1e-7);
214 }
215 ASSERT_LE(d, 1e-7);
216
217 spqlios_free(a2);
218 spqlios_free(b2);
219 // spqlios_free(c2);
220}
221
222TEST(fft, fftvec_convolution_recursiveoverk) {
223 static const uint64_t nn = 32768; // vary accross (8192, 16384), 32768, 65536
224 double* a = (double*)spqlios_alloc_custom_align(32, nn * 8);
225 double* a2 = (double*)spqlios_alloc_custom_align(32, nn * 8);
226 double* b = (double*)spqlios_alloc_custom_align(32, nn * 8);
227 double* dist_vector = (double*)spqlios_alloc_custom_align(32, nn * 8);
228
229 printf("N size: %" PRId64 "\n", nn);
230
231 for (uint32_t k = 14; k <= 24; k++) { // vary k
232 printf("k size: %" PRId32 "\n", k);
233 int64_t p = UINT64_C(1) << k;
234 for (uint32_t i = 0; i < nn; i++) {
235 a[i] = (rand() % p) - p / 2;
236 b[i] = (rand() % p) - p / 2;
237 a2[i] = 0;
238 }
239 cplx_fft_simple(nn / 2, a);
240 cplx_fft_simple(nn / 2, b);
241 cplx_fftvec_addmul_simple(nn / 2, a2, a, b);
242 cplx_ifft_simple(nn / 2, a2);
243 double distance = 0;
244 for (uint32_t i = 0; i < nn; i++) {
245 double curdist = fabs(a2[i] / (nn / 2.) - rint(a2[i] / (nn / 2.)));
246 if (distance < curdist) distance = curdist;
247 dist_vector[i] = a2[i] / (nn / 2.) - rint(a2[i] / (nn / 2.));
248 }
249 printf("distance: %lf\n", distance);
250 ASSERT_LE(distance, 0.5); // switch from previous 0.1 to 0.5 per experiment 1 reqs
251 double mean = 0;
252 for (uint32_t i = 0; i < nn; i++) {
253 mean = mean + dist_vector[i];
254 }
255 mean = mean / nn;
256 double variance = 0;
257 for (uint32_t i = 0; i < nn; i++) {
258 variance = variance + pow((mean - dist_vector[i]), 2);
259 }
260 double stdev = sqrt(variance / nn);
261 printf("stdev: %lf\n", stdev);
262 }
263
264 spqlios_free(a);
265 spqlios_free(b);
266 spqlios_free(a2);
267 spqlios_free(dist_vector);
268}
269
270#ifdef __x86_64__
271TEST(fft, cplx_fft_ref_vs_fft_reim_ref) {
272 for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
273 uint64_t m = nn / 2;
274 CPLX_FFT_PRECOMP* tables = new_cplx_fft_precomp(m, 0);
275 REIM_FFT_PRECOMP* reimtables = new_reim_fft_precomp(m, 0);
276 CPLX* a = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
277 CPLX* a1 = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
278 double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
279 int64_t p = 1 << 16;
280 for (uint32_t i = 0; i < nn / 2; i++) {
281 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
282 a[i][1] = (rand() % p) - p / 2;
283 }
284 memcpy(a1, a, nn / 2 * sizeof(CPLX));
285 for (uint32_t i = 0; i < nn / 2; i++) {
286 a2[i] = a[i][0];
287 a2[nn / 2 + i] = a[i][1];
288 }
289 cplx_fft_ref(tables, a1);
290 reim_fft_ref(reimtables, a2);
291 double d = 0;
292 for (uint32_t i = 0; i < nn / 2; i++) {
293 double dre = fabs(a1[i][0] - a2[i]);
294 double dim = fabs(a1[i][1] - a2[nn / 2 + i]);
295 if (dre > d) d = dre;
296 if (dim > d) d = dim;
297 ASSERT_LE(d, 1e-7);
298 }
299 ASSERT_LE(d, 1e-7);
300 spqlios_free(a);
301 spqlios_free(a1);
302 spqlios_free(a2);
303 delete_cplx_fft_precomp(tables);
304 delete_reim_fft_precomp(reimtables);
305 }
306}
307#endif
308
309TEST(fft, cplx_ifft_ref_vs_reim_ifft_ref) {
310 for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
311 uint64_t m = nn / 2;
312 CPLX_IFFT_PRECOMP* tables = new_cplx_ifft_precomp(m, 0);
313 REIM_IFFT_PRECOMP* reimtables = new_reim_ifft_precomp(m, 0);
314 CPLX* a = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
315 CPLX* a1 = (CPLX*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
316 double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
317 int64_t p = 1 << 16;
318 for (uint32_t i = 0; i < nn / 2; i++) {
319 a[i][0] = (rand() % p) - p / 2; // between -p/2 and p/2
320 a[i][1] = (rand() % p) - p / 2;
321 }
322 memcpy(a1, a, nn / 2 * sizeof(CPLX));
323 for (uint32_t i = 0; i < nn / 2; i++) {
324 a2[i] = a[i][0];
325 a2[nn / 2 + i] = a[i][1];
326 }
327 cplx_ifft_ref(tables, a1);
328 reim_ifft_ref(reimtables, a2);
329 double d = 0;
330 for (uint32_t i = 0; i < nn / 2; i++) {
331 double dre = fabs(a1[i][0] - a2[i]);
332 double dim = fabs(a1[i][1] - a2[nn / 2 + i]);
333 if (dre > d) d = dre;
334 if (dim > d) d = dim;
335 ASSERT_LE(d, 1e-7);
336 }
337 ASSERT_LE(d, 1e-7);
338 spqlios_free(a);
339 spqlios_free(a1);
340 spqlios_free(a2);
341 delete_cplx_fft_precomp(tables);
342 delete_reim_fft_precomp(reimtables);
343 }
344}
345
346#ifdef __x86_64__
347TEST(fft, reim4_vecfft_addmul_fma_vs_ref) {
348 for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
349 uint64_t m = nn / 2;
350 REIM4_FFTVEC_ADDMUL_PRECOMP* tbl = new_reim4_fftvec_addmul_precomp(m);
351 ASSERT_TRUE(tbl != nullptr);
352 double* a1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
353 double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
354 double* b1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
355 double* b2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
356 double* r1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
357 double* r2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
358 int64_t p = 1 << 16;
359 for (uint32_t i = 0; i < nn; i++) {
360 a1[i] = (rand() % p) - p / 2; // between -p/2 and p/2
361 b1[i] = (rand() % p) - p / 2;
362 r1[i] = (rand() % p) - p / 2;
363 }
364 memcpy(a2, a1, nn / 2 * sizeof(CPLX));
365 memcpy(b2, b1, nn / 2 * sizeof(CPLX));
366 memcpy(r2, r1, nn / 2 * sizeof(CPLX));
367 reim4_fftvec_addmul_ref(tbl, r1, a1, b1);
368 reim4_fftvec_addmul_fma(tbl, r2, a2, b2);
369 double d = 0;
370 for (uint32_t i = 0; i < nn; i++) {
371 double di = fabs(r1[i] - r2[i]);
372 if (di > d) d = di;
373 ASSERT_LE(d, 1e-8);
374 }
375 ASSERT_LE(d, 1e-8);
376 spqlios_free(a1);
377 spqlios_free(a2);
378 spqlios_free(b1);
379 spqlios_free(b2);
380 spqlios_free(r1);
381 spqlios_free(r2);
382 delete_reim4_fftvec_addmul_precomp(tbl);
383 }
384}
385#endif
386
387#ifdef __x86_64__
388TEST(fft, reim4_vecfft_mul_fma_vs_ref) {
389 for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
390 uint64_t m = nn / 2;
391 REIM4_FFTVEC_MUL_PRECOMP* tbl = new_reim4_fftvec_mul_precomp(m);
392 double* a1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
393 double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
394 double* b1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
395 double* b2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
396 double* r1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
397 double* r2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
398 int64_t p = 1 << 16;
399 for (uint32_t i = 0; i < nn; i++) {
400 a1[i] = (rand() % p) - p / 2; // between -p/2 and p/2
401 b1[i] = (rand() % p) - p / 2;
402 r1[i] = (rand() % p) - p / 2;
403 }
404 memcpy(a2, a1, nn / 2 * sizeof(CPLX));
405 memcpy(b2, b1, nn / 2 * sizeof(CPLX));
406 memcpy(r2, r1, nn / 2 * sizeof(CPLX));
407 reim4_fftvec_mul_ref(tbl, r1, a1, b1);
408 reim4_fftvec_mul_fma(tbl, r2, a2, b2);
409 double d = 0;
410 for (uint32_t i = 0; i < nn; i++) {
411 double di = fabs(r1[i] - r2[i]);
412 if (di > d) d = di;
413 ASSERT_LE(d, 1e-8);
414 }
415 ASSERT_LE(d, 1e-8);
416 spqlios_free(a1);
417 spqlios_free(a2);
418 spqlios_free(b1);
419 spqlios_free(b2);
420 spqlios_free(r1);
421 spqlios_free(r2);
422 delete_reim_fftvec_mul_precomp(tbl);
423 }
424}
425#endif
426
427#ifdef __x86_64__
428TEST(fft, reim4_from_cplx_fma_vs_ref) {
429 for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
430 uint64_t m = nn / 2;
431 REIM4_FROM_CPLX_PRECOMP* tbl = new_reim4_from_cplx_precomp(m);
432 double* a1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
433 double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
434 double* r1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
435 double* r2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
436 int64_t p = 1 << 16;
437 for (uint32_t i = 0; i < nn; i++) {
438 a1[i] = (rand() % p) - p / 2; // between -p/2 and p/2
439 r1[i] = (rand() % p) - p / 2;
440 }
441 memcpy(a2, a1, nn / 2 * sizeof(CPLX));
442 memcpy(r2, r1, nn / 2 * sizeof(CPLX));
443 reim4_from_cplx_ref(tbl, r1, a1);
444 reim4_from_cplx_fma(tbl, r2, a2);
445 double d = 0;
446 for (uint32_t i = 0; i < nn; i++) {
447 double di = fabs(r1[i] - r2[i]);
448 if (di > d) d = di;
449 ASSERT_LE(d, 1e-8);
450 }
451 ASSERT_LE(d, 1e-8);
452 spqlios_free(a1);
453 spqlios_free(a2);
454 spqlios_free(r1);
455 spqlios_free(r2);
456 delete_reim4_from_cplx_precomp(tbl);
457 }
458}
459#endif
460
461#ifdef __x86_64__
462TEST(fft, reim4_to_cplx_fma_vs_ref) {
463 for (uint64_t nn : {16, 32, 64, 1024, 8192, 65536}) {
464 uint64_t m = nn / 2;
465 REIM4_TO_CPLX_PRECOMP* tbl = new_reim4_to_cplx_precomp(m);
466 double* a1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
467 double* a2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
468 double* r1 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
469 double* r2 = (double*)spqlios_alloc_custom_align(32, nn / 2 * sizeof(CPLX));
470 int64_t p = 1 << 16;
471 for (uint32_t i = 0; i < nn; i++) {
472 a1[i] = (rand() % p) - p / 2; // between -p/2 and p/2
473 r1[i] = (rand() % p) - p / 2;
474 }
475 memcpy(a2, a1, nn / 2 * sizeof(CPLX));
476 memcpy(r2, r1, nn / 2 * sizeof(CPLX));
477 reim4_to_cplx_ref(tbl, r1, a1);
478 reim4_to_cplx_fma(tbl, r2, a2);
479 double d = 0;
480 for (uint32_t i = 0; i < nn; i++) {
481 double di = fabs(r1[i] - r2[i]);
482 if (di > d) d = di;
483 ASSERT_LE(d, 1e-8);
484 }
485 ASSERT_LE(d, 1e-8);
486 spqlios_free(a1);
487 spqlios_free(a2);
488 spqlios_free(r1);
489 spqlios_free(r2);
490 delete_reim4_from_cplx_precomp(tbl);
491 }
492}
493#endif