poulpy/poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/reim/reim_ifft_ref.c
poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/reim/reim_ifft_ref.cRaw
1#include "../commons_private.h"
2#include "reim_fft_internal.h"
3#include "reim_fft_private.h"
4
5void reim_invctwiddle(double* ra, double* ia, double* rb, double* ib, double omre, double omim) {
6 double rdiff = *ra - *rb;
7 double idiff = *ia - *ib;
8 *ra = *ra + *rb;
9 *ia = *ia + *ib;
10 *rb = rdiff * omre - idiff * omim;
11 *ib = rdiff * omim + idiff * omre;
12}
13
14// i (omre + i omim) = -omim + i omre
15void reim_invcitwiddle(double* ra, double* ia, double* rb, double* ib, double omre, double omim) {
16 double rdiff = *ra - *rb;
17 double idiff = *ia - *ib;
18 *ra = *ra + *rb;
19 *ia = *ia + *ib;
20 *rb = rdiff * omim + idiff * omre;
21 *ib = -rdiff * omre + idiff * omim;
22}
23
24void reim_ifft16_ref(double* dre, double* dim, const void* pom) {
25 const double* om = (const double*)pom;
26 {
27 double omare = om[0];
28 double ombre = om[1];
29 double omcre = om[2];
30 double omdre = om[3];
31 double omaim = om[4];
32 double ombim = om[5];
33 double omcim = om[6];
34 double omdim = om[7];
35 reim_invctwiddle(&dre[0], &dim[0], &dre[1], &dim[1], omare, omaim);
36 reim_invcitwiddle(&dre[2], &dim[2], &dre[3], &dim[3], omare, omaim);
37 reim_invctwiddle(&dre[4], &dim[4], &dre[5], &dim[5], ombre, ombim);
38 reim_invcitwiddle(&dre[6], &dim[6], &dre[7], &dim[7], ombre, ombim);
39 reim_invctwiddle(&dre[8], &dim[8], &dre[9], &dim[9], omcre, omcim);
40 reim_invcitwiddle(&dre[10], &dim[10], &dre[11], &dim[11], omcre, omcim);
41 reim_invctwiddle(&dre[12], &dim[12], &dre[13], &dim[13], omdre, omdim);
42 reim_invcitwiddle(&dre[14], &dim[14], &dre[15], &dim[15], omdre, omdim);
43 }
44 {
45 double omare = om[8];
46 double omaim = om[9];
47 double ombre = om[10];
48 double ombim = om[11];
49 reim_invctwiddle(&dre[0], &dim[0], &dre[2], &dim[2], omare, omaim);
50 reim_invctwiddle(&dre[1], &dim[1], &dre[3], &dim[3], omare, omaim);
51 reim_invcitwiddle(&dre[4], &dim[4], &dre[6], &dim[6], omare, omaim);
52 reim_invcitwiddle(&dre[5], &dim[5], &dre[7], &dim[7], omare, omaim);
53 reim_invctwiddle(&dre[8], &dim[8], &dre[10], &dim[10], ombre, ombim);
54 reim_invctwiddle(&dre[9], &dim[9], &dre[11], &dim[11], ombre, ombim);
55 reim_invcitwiddle(&dre[12], &dim[12], &dre[14], &dim[14], ombre, ombim);
56 reim_invcitwiddle(&dre[13], &dim[13], &dre[15], &dim[15], ombre, ombim);
57 }
58 {
59 double omre = om[12];
60 double omim = om[13];
61 reim_invctwiddle(&dre[0], &dim[0], &dre[4], &dim[4], omre, omim);
62 reim_invctwiddle(&dre[1], &dim[1], &dre[5], &dim[5], omre, omim);
63 reim_invctwiddle(&dre[2], &dim[2], &dre[6], &dim[6], omre, omim);
64 reim_invctwiddle(&dre[3], &dim[3], &dre[7], &dim[7], omre, omim);
65 reim_invcitwiddle(&dre[8], &dim[8], &dre[12], &dim[12], omre, omim);
66 reim_invcitwiddle(&dre[9], &dim[9], &dre[13], &dim[13], omre, omim);
67 reim_invcitwiddle(&dre[10], &dim[10], &dre[14], &dim[14], omre, omim);
68 reim_invcitwiddle(&dre[11], &dim[11], &dre[15], &dim[15], omre, omim);
69 }
70 {
71 double omre = om[14];
72 double omim = om[15];
73 reim_invctwiddle(&dre[0], &dim[0], &dre[8], &dim[8], omre, omim);
74 reim_invctwiddle(&dre[1], &dim[1], &dre[9], &dim[9], omre, omim);
75 reim_invctwiddle(&dre[2], &dim[2], &dre[10], &dim[10], omre, omim);
76 reim_invctwiddle(&dre[3], &dim[3], &dre[11], &dim[11], omre, omim);
77 reim_invctwiddle(&dre[4], &dim[4], &dre[12], &dim[12], omre, omim);
78 reim_invctwiddle(&dre[5], &dim[5], &dre[13], &dim[13], omre, omim);
79 reim_invctwiddle(&dre[6], &dim[6], &dre[14], &dim[14], omre, omim);
80 reim_invctwiddle(&dre[7], &dim[7], &dre[15], &dim[15], omre, omim);
81 }
82}
83
84void fill_reim_ifft16_omegas(const double entry_pwr, double** omg) {
85 const double j_pow = 1. / 8.;
86 const double k_pow = 1. / 16.;
87 const double pin = entry_pwr / 2.;
88 const double pin_2 = entry_pwr / 4.;
89 const double pin_4 = entry_pwr / 8.;
90 const double pin_8 = entry_pwr / 16.;
91 // ((8,9,10,11),(12,13,14,15)) are 4 reals then 4 imag of om^1/8*(1,k,j,kj)
92 (*omg)[0] = cos(2. * M_PI * (pin_8));
93 (*omg)[1] = cos(2. * M_PI * (pin_8 + j_pow));
94 (*omg)[2] = cos(2. * M_PI * (pin_8 + k_pow));
95 (*omg)[3] = cos(2. * M_PI * (pin_8 + j_pow + k_pow));
96 (*omg)[4] = -sin(2. * M_PI * (pin_8));
97 (*omg)[5] = -sin(2. * M_PI * (pin_8 + j_pow));
98 (*omg)[6] = -sin(2. * M_PI * (pin_8 + k_pow));
99 (*omg)[7] = -sin(2. * M_PI * (pin_8 + j_pow + k_pow));
100 // (4,5) and (6,7) are real and imag of om^1/4 and j.om^1/4
101 (*omg)[8] = cos(2. * M_PI * (pin_4));
102 (*omg)[9] = -sin(2. * M_PI * (pin_4));
103 (*omg)[10] = cos(2. * M_PI * (pin_4 + j_pow));
104 (*omg)[11] = -sin(2. * M_PI * (pin_4 + j_pow));
105 // 2 and 3 are real and imag of om^1/2
106 (*omg)[12] = cos(2. * M_PI * (pin_2));
107 (*omg)[13] = -sin(2. * M_PI * (pin_2));
108 // 0 and 1 are real and imag of om
109 (*omg)[14] = cos(2. * M_PI * pin);
110 (*omg)[15] = -sin(2. * M_PI * pin);
111 *omg += 16;
112}
113
114void reim_ifft8_ref(double* dre, double* dim, const void* pom) {
115 const double* om = (const double*)pom;
116 {
117 double omare = om[0];
118 double ombre = om[1];
119 double omaim = om[2];
120 double ombim = om[3];
121 reim_invctwiddle(&dre[0], &dim[0], &dre[1], &dim[1], omare, omaim);
122 reim_invcitwiddle(&dre[2], &dim[2], &dre[3], &dim[3], omare, omaim);
123 reim_invctwiddle(&dre[4], &dim[4], &dre[5], &dim[5], ombre, ombim);
124 reim_invcitwiddle(&dre[6], &dim[6], &dre[7], &dim[7], ombre, ombim);
125 }
126 {
127 double omare = om[4];
128 double omaim = om[5];
129 reim_invctwiddle(&dre[0], &dim[0], &dre[2], &dim[2], omare, omaim);
130 reim_invctwiddle(&dre[1], &dim[1], &dre[3], &dim[3], omare, omaim);
131 reim_invcitwiddle(&dre[4], &dim[4], &dre[6], &dim[6], omare, omaim);
132 reim_invcitwiddle(&dre[5], &dim[5], &dre[7], &dim[7], omare, omaim);
133 }
134 {
135 double omre = om[6];
136 double omim = om[7];
137 reim_invctwiddle(&dre[0], &dim[0], &dre[4], &dim[4], omre, omim);
138 reim_invctwiddle(&dre[1], &dim[1], &dre[5], &dim[5], omre, omim);
139 reim_invctwiddle(&dre[2], &dim[2], &dre[6], &dim[6], omre, omim);
140 reim_invctwiddle(&dre[3], &dim[3], &dre[7], &dim[7], omre, omim);
141 }
142}
143
144void fill_reim_ifft8_omegas(const double entry_pwr, double** omg) {
145 const double j_pow = 1. / 8.;
146 const double pin = entry_pwr / 2.;
147 const double pin_2 = entry_pwr / 4.;
148 const double pin_4 = entry_pwr / 8.;
149 // (4,5) and (6,7) are real and imag of om^1/4 and j.om^1/4
150 (*omg)[0] = cos(2. * M_PI * (pin_4));
151 (*omg)[1] = cos(2. * M_PI * (pin_4 + j_pow));
152 (*omg)[2] = -sin(2. * M_PI * (pin_4));
153 (*omg)[3] = -sin(2. * M_PI * (pin_4 + j_pow));
154 // 2 and 3 are real and imag of om^1/2
155 (*omg)[4] = cos(2. * M_PI * (pin_2));
156 (*omg)[5] = -sin(2. * M_PI * (pin_2));
157 // 0 and 1 are real and imag of om
158 (*omg)[6] = cos(2. * M_PI * pin);
159 (*omg)[7] = -sin(2. * M_PI * pin);
160 *omg += 8;
161}
162
163void reim_ifft4_ref(double* dre, double* dim, const void* pom) {
164 const double* om = (const double*)pom;
165 {
166 double omare = om[0];
167 double omaim = om[1];
168 reim_invctwiddle(&dre[0], &dim[0], &dre[1], &dim[1], omare, omaim);
169 reim_invcitwiddle(&dre[2], &dim[2], &dre[3], &dim[3], omare, omaim);
170 }
171 {
172 double omare = om[2];
173 double omaim = om[3];
174 reim_invctwiddle(&dre[0], &dim[0], &dre[2], &dim[2], omare, omaim);
175 reim_invctwiddle(&dre[1], &dim[1], &dre[3], &dim[3], omare, omaim);
176 }
177}
178
179void fill_reim_ifft4_omegas(const double entry_pwr, double** omg) {
180 const double pin = entry_pwr / 2.;
181 const double pin_2 = entry_pwr / 4.;
182 // 2 and 3 are real and imag of om^1/2
183 (*omg)[0] = cos(2. * M_PI * (pin_2));
184 (*omg)[1] = -sin(2. * M_PI * (pin_2));
185 // 0 and 1 are real and imag of om
186 (*omg)[2] = cos(2. * M_PI * pin);
187 (*omg)[3] = -sin(2. * M_PI * pin);
188 *omg += 4;
189}
190
191void reim_ifft2_ref(double* dre, double* dim, const void* pom) {
192 const double* om = (const double*)pom;
193 {
194 double omare = om[0];
195 double omaim = om[1];
196 reim_invctwiddle(&dre[0], &dim[0], &dre[1], &dim[1], omare, omaim);
197 }
198}
199
200void fill_reim_ifft2_omegas(const double entry_pwr, double** omg) {
201 const double pin = entry_pwr / 2.;
202 // 0 and 1 are real and imag of om
203 (*omg)[0] = cos(2. * M_PI * pin);
204 (*omg)[1] = -sin(2. * M_PI * pin);
205 *omg += 2;
206}
207
208void reim_invtwiddle_ifft_ref(uint64_t h, double* re, double* im, double om[2]) {
209 for (uint64_t i = 0; i < h; ++i) {
210 reim_invctwiddle(&re[i], &im[i], &re[h + i], &im[h + i], om[0], om[1]);
211 }
212}
213
214void reim_invbitwiddle_ifft_ref(uint64_t h, double* re, double* im, double om[4]) {
215 double* r0 = re;
216 double* r1 = re + h;
217 double* r2 = re + 2 * h;
218 double* r3 = re + 3 * h;
219 double* i0 = im;
220 double* i1 = im + h;
221 double* i2 = im + 2 * h;
222 double* i3 = im + 3 * h;
223 for (uint64_t i = 0; i < h; ++i) {
224 reim_invctwiddle(&r0[i], &i0[i], &r1[i], &i1[i], om[0], om[1]);
225 reim_invcitwiddle(&r2[i], &i2[i], &r3[i], &i3[i], om[0], om[1]);
226 }
227 for (uint64_t i = 0; i < h; ++i) {
228 reim_invctwiddle(&r0[i], &i0[i], &r2[i], &i2[i], om[2], om[3]);
229 reim_invctwiddle(&r1[i], &i1[i], &r3[i], &i3[i], om[2], om[3]);
230 }
231}
232
233void reim_ifft_bfs_16_ref(uint64_t m, double* re, double* im, double** omg) {
234 uint64_t log2m = log2(m);
235 for (uint64_t off = 0; off < m; off += 16) {
236 reim_ifft16_ref(re + off, im + off, *omg);
237 *omg += 16;
238 }
239 uint64_t h = 16;
240 uint64_t ms2 = m / 2;
241 while (h < ms2) {
242 uint64_t mm = h << 2;
243 for (uint64_t off = 0; off < m; off += mm) {
244 reim_invbitwiddle_ifft_ref(h, re + off, im + off, *omg);
245 *omg += 4;
246 }
247 h = mm;
248 }
249 if (log2m % 2 != 0) {
250 if (h != ms2) abort(); // bug
251 // do the first twiddle iteration normally
252 reim_invtwiddle_ifft_ref(h, re, im, *omg);
253 *omg += 2;
254 h = m;
255 }
256}
257
258void fill_reim_ifft_bfs_16_omegas(uint64_t m, double entry_pwr, double** omg) {
259 uint64_t log2m = log2(m);
260 // uint64_t mm = 16;
261 double ss = entry_pwr * 16. / m;
262 for (uint64_t off = 0; off < m; off += 16) {
263 double s = ss + fracrevbits(off / 16);
264 fill_reim_ifft16_omegas(s, omg);
265 }
266 uint64_t h = 16;
267 uint64_t ms2 = m / 2;
268 while (h < ms2) {
269 uint64_t mm = h << 2;
270 for (uint64_t off = 0; off < m; off += mm) {
271 double rs0 = ss + fracrevbits(off / mm) / 4.;
272 double rs1 = 2. * rs0;
273 (*omg)[0] = cos(2 * M_PI * rs0);
274 (*omg)[1] = -sin(2 * M_PI * rs0);
275 (*omg)[2] = cos(2 * M_PI * rs1);
276 (*omg)[3] = -sin(2 * M_PI * rs1);
277 *omg += 4;
278 }
279 ss *= 4.;
280 h = mm;
281 }
282 if (log2m % 2 != 0) {
283 if (h != ms2) abort(); // bug
284 // do the first twiddle iteration normally
285 (*omg)[0] = cos(2 * M_PI * ss);
286 (*omg)[1] = -sin(2 * M_PI * ss);
287 *omg += 2;
288 h = m;
289 ss *= 2.;
290 }
291 if (ss != entry_pwr) abort();
292}
293
294void reim_ifft_rec_16_ref(uint64_t m, double* re, double* im, double** omg) {
295 if (m <= 2048) return reim_ifft_bfs_16_ref(m, re, im, omg);
296 const uint32_t h = m / 2;
297 reim_ifft_rec_16_ref(h, re, im, omg);
298 reim_ifft_rec_16_ref(h, re + h, im + h, omg);
299 reim_invtwiddle_ifft_ref(h, re, im, *omg);
300 *omg += 2;
301}
302
303void fill_reim_ifft_rec_16_omegas(uint64_t m, double entry_pwr, double** omg) {
304 if (m <= 2048) return fill_reim_ifft_bfs_16_omegas(m, entry_pwr, omg);
305 const uint64_t h = m / 2;
306 const double s = entry_pwr / 2;
307 fill_reim_ifft_rec_16_omegas(h, s, omg);
308 fill_reim_ifft_rec_16_omegas(h, s + 0.5, omg);
309 (*omg)[0] = cos(2 * M_PI * s);
310 (*omg)[1] = -sin(2 * M_PI * s);
311 *omg += 2;
312}
313
314void reim_ifft_ref(const REIM_IFFT_PRECOMP* precomp, double* dat) {
315 const int32_t m = precomp->m;
316 double* omg = precomp->powomegas;
317 double* re = dat;
318 double* im = dat + m;
319 if (m <= 16) {
320 switch (m) {
321 case 1:
322 return;
323 case 2:
324 return reim_ifft2_ref(re, im, omg);
325 case 4:
326 return reim_ifft4_ref(re, im, omg);
327 case 8:
328 return reim_ifft8_ref(re, im, omg);
329 case 16:
330 return reim_ifft16_ref(re, im, omg);
331 default:
332 abort(); // m is not a power of 2
333 }
334 }
335 if (m <= 2048) return reim_ifft_bfs_16_ref(m, re, im, &omg);
336 return reim_ifft_rec_16_ref(m, re, im, &omg);
337}
338
339EXPORT REIM_IFFT_PRECOMP* new_reim_ifft_precomp(uint32_t m, uint32_t num_buffers) {
340 const uint64_t OMG_SPACE = ceilto64b(2 * m * sizeof(double));
341 const uint64_t BUF_SIZE = ceilto64b(2 * m * sizeof(double));
342 void* reps = malloc(sizeof(REIM_IFFT_PRECOMP) // base
343 + 63 // padding
344 + OMG_SPACE // tables //TODO 16?
345 + num_buffers * BUF_SIZE // buffers
346 );
347 uint64_t aligned_addr = ceilto64b((uint64_t)(reps) + sizeof(REIM_IFFT_PRECOMP));
348 REIM_IFFT_PRECOMP* r = (REIM_IFFT_PRECOMP*)reps;
349 r->m = m;
350 r->buf_size = BUF_SIZE;
351 r->powomegas = (double*)aligned_addr;
352 r->aligned_buffers = (void*)(aligned_addr + OMG_SPACE);
353 // fill in powomegas
354 double* omg = (double*)r->powomegas;
355 if (m <= 16) {
356 switch (m) {
357 case 1:
358 break;
359 case 2:
360 fill_reim_ifft2_omegas(0.25, &omg);
361 break;
362 case 4:
363 fill_reim_ifft4_omegas(0.25, &omg);
364 break;
365 case 8:
366 fill_reim_ifft8_omegas(0.25, &omg);
367 break;
368 case 16:
369 fill_reim_ifft16_omegas(0.25, &omg);
370 break;
371 default:
372 abort(); // m is not a power of 2
373 }
374 } else if (m <= 2048) {
375 fill_reim_ifft_bfs_16_omegas(m, 0.25, &omg);
376 } else {
377 fill_reim_ifft_rec_16_omegas(m, 0.25, &omg);
378 }
379 if (((uint64_t)omg) - aligned_addr > OMG_SPACE) abort();
380 // dispatch the right implementation
381 {
382 if (CPU_SUPPORTS("fma")) {
383 r->function = reim_ifft_avx2_fma;
384 } else {
385 r->function = reim_ifft_ref;
386 }
387 }
388 return reps;
389}
390
391void reim_naive_ifft(uint64_t m, double entry_pwr, double* re, double* im) {
392 if (m == 1) return;
393 // twiddle
394 const uint64_t h = m / 2;
395 const double s = entry_pwr / 2.;
396 reim_naive_ifft(h, s, re, im);
397 reim_naive_ifft(h, s + 0.5, re + h, im + h);
398 const double sre = cos(2 * M_PI * s);
399 const double sim = -sin(2 * M_PI * s);
400 for (uint64_t j = 0; j < h; ++j) {
401 double rdiff = re[j] - re[h + j];
402 double idiff = im[j] - im[h + j];
403 re[j] = re[j] + re[h + j];
404 im[j] = im[j] + im[h + j];
405 re[h + j] = rdiff * sre - idiff * sim;
406 im[h + j] = idiff * sre + rdiff * sim;
407 }
408}