poulpy/poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/reim/reim_fft_neon.c
poulpy-backend/src/cpu_spqlios/spqlios-arithmetic/spqlios/reim/reim_fft_neon.cRaw
1/*
2 * This file is adapted from the implementation of the FFT on Arm64/Neon
3 * available in https://github.com/cothan/Falcon-Arm (neon/fft.c).
4 * =============================================================================
5 * Copyright (c) 2022 by Cryptographic Engineering Research Group (CERG)
6 * ECE Department, George Mason University
7 * Fairfax, VA, U.S.A.
8 * @author: Duc Tri Nguyen dnguye69@gmu.edu, cothannguyen@gmail.com
9 * Licensed under the Apache License, Version 2.0 (the "License");
10 * =============================================================================
11 *
12 * The original source file has been modified by the authors of spqlios-arithmetic
13 * to be interfaced with the dynamic twiddle factors generator of the spqlios-arithmetic
14 * library, as well as the recursive dfs strategy for large complex dimensions m>=2048.
15 */
16
17#include <math.h>
18
19#include "../commons.h"
20#include "../ext/neon_accel/macrof.h"
21#include "../ext/neon_accel/macrofx4.h"
22
23void fill_reim_fft16_omegas_neon(const double entry_pwr, double** omg) {
24 const double j_pow = 1. / 8.;
25 const double k_pow = 1. / 16.;
26 const double pin = entry_pwr / 2.;
27 const double pin_2 = entry_pwr / 4.;
28 const double pin_4 = entry_pwr / 8.;
29 const double pin_8 = entry_pwr / 16.;
30 // 0 and 1 are real and imag of om
31 (*omg)[0] = cos(2. * M_PI * pin);
32 (*omg)[1] = sin(2. * M_PI * pin);
33 // 2 and 3 are real and imag of om^1/2
34 (*omg)[2] = cos(2. * M_PI * (pin_2));
35 (*omg)[3] = sin(2. * M_PI * (pin_2));
36 // (4,5) and (6,7) are real and imag of om^1/4 and j.om^1/4
37 (*omg)[4] = cos(2. * M_PI * (pin_4));
38 (*omg)[5] = sin(2. * M_PI * (pin_4));
39 (*omg)[6] = cos(2. * M_PI * (pin_4 + j_pow));
40 (*omg)[7] = sin(2. * M_PI * (pin_4 + j_pow));
41 // ((8,9,10,11),(12,13,14,15)) are 4 reals then 4 imag of om^1/8*(1,k,j,kj)
42 (*omg)[8] = cos(2. * M_PI * (pin_8));
43 (*omg)[10] = cos(2. * M_PI * (pin_8 + j_pow));
44 (*omg)[12] = cos(2. * M_PI * (pin_8 + k_pow));
45 (*omg)[14] = cos(2. * M_PI * (pin_8 + j_pow + k_pow));
46 (*omg)[9] = sin(2. * M_PI * (pin_8));
47 (*omg)[11] = sin(2. * M_PI * (pin_8 + j_pow));
48 (*omg)[13] = sin(2. * M_PI * (pin_8 + k_pow));
49 (*omg)[15] = sin(2. * M_PI * (pin_8 + j_pow + k_pow));
50 *omg += 16;
51}
52
53EXPORT void reim_fft16_neon(double* dre, double* dim, const void* omega) {
54 const double* pom = omega;
55 // Total SIMD register: 28 = 24 + 4
56 float64x2x2_t s_re_im; // 2
57 float64x2x4_t x_re, x_im, y_re, y_im, t_re, t_im, v_re, v_im; // 32
58
59 {
60 /*
61 Level 2
62 ( 8, 24) * ( 0, 1)
63 ( 9, 25) * ( 0, 1)
64 ( 10, 26) * ( 0, 1)
65 ( 11, 27) * ( 0, 1)
66 ( 12, 28) * ( 0, 1)
67 ( 13, 29) * ( 0, 1)
68 ( 14, 30) * ( 0, 1)
69 ( 15, 31) * ( 0, 1)
70
71 ( 8, 24) = ( 0, 16) - @
72 ( 9, 25) = ( 1, 17) - @
73 ( 10, 26) = ( 2, 18) - @
74 ( 11, 27) = ( 3, 19) - @
75 ( 12, 28) = ( 4, 20) - @
76 ( 13, 29) = ( 5, 21) - @
77 ( 14, 30) = ( 6, 22) - @
78 ( 15, 31) = ( 7, 23) - @
79
80 ( 0, 16) = ( 0, 16) + @
81 ( 1, 17) = ( 1, 17) + @
82 ( 2, 18) = ( 2, 18) + @
83 ( 3, 19) = ( 3, 19) + @
84 ( 4, 20) = ( 4, 20) + @
85 ( 5, 21) = ( 5, 21) + @
86 ( 6, 22) = ( 6, 22) + @
87 ( 7, 23) = ( 7, 23) + @
88 */
89 vload(s_re_im.val[0], pom);
90
91 vloadx4(y_re, dre + 8);
92 vloadx4(y_im, dim + 8);
93
94 FWD_TOP_LANEx4(v_re, v_im, y_re, y_im, s_re_im.val[0]);
95
96 vloadx4(x_re, dre);
97 vloadx4(x_im, dim);
98
99 FWD_BOTx4(x_re, x_im, y_re, y_im, v_re, v_im);
100
101 // vstorex4(dre, x_re);
102 // vstorex4(dim, x_im);
103 // vstorex4(dre + 8, y_re);
104 // vstorex4(dim + 8, y_im);
105 // return;
106 }
107 {
108 /*
109 Level 3
110
111 ( 4, 20) * ( 0, 1)
112 ( 5, 21) * ( 0, 1)
113 ( 6, 22) * ( 0, 1)
114 ( 7, 23) * ( 0, 1)
115
116 ( 4, 20) = ( 0, 16) - @
117 ( 5, 21) = ( 1, 17) - @
118 ( 6, 22) = ( 2, 18) - @
119 ( 7, 23) = ( 3, 19) - @
120
121 ( 0, 16) = ( 0, 16) + @
122 ( 1, 17) = ( 1, 17) + @
123 ( 2, 18) = ( 2, 18) + @
124 ( 3, 19) = ( 3, 19) + @
125
126 ( 12, 28) * ( 0, 1)
127 ( 13, 29) * ( 0, 1)
128 ( 14, 30) * ( 0, 1)
129 ( 15, 31) * ( 0, 1)
130
131 ( 12, 28) = ( 8, 24) - j@
132 ( 13, 29) = ( 9, 25) - j@
133 ( 14, 30) = ( 10, 26) - j@
134 ( 15, 31) = ( 11, 27) - j@
135
136 ( 8, 24) = ( 8, 24) + j@
137 ( 9, 25) = ( 9, 25) + j@
138 ( 10, 26) = ( 10, 26) + j@
139 ( 11, 27) = ( 11, 27) + j@
140 */
141
142 // vloadx4(y_re, dre + 8);
143 // vloadx4(y_im, dim + 8);
144 // vloadx4(x_re, dre);
145 // vloadx4(x_im, dim);
146
147 vload(s_re_im.val[0], pom + 2);
148
149 FWD_TOP_LANE(t_re.val[0], t_im.val[0], x_re.val[2], x_im.val[2], s_re_im.val[0]);
150 FWD_TOP_LANE(t_re.val[1], t_im.val[1], x_re.val[3], x_im.val[3], s_re_im.val[0]);
151 FWD_TOP_LANE(t_re.val[2], t_im.val[2], y_re.val[2], y_im.val[2], s_re_im.val[0]);
152 FWD_TOP_LANE(t_re.val[3], t_im.val[3], y_re.val[3], y_im.val[3], s_re_im.val[0]);
153
154 FWD_BOT(x_re.val[0], x_im.val[0], x_re.val[2], x_im.val[2], t_re.val[0], t_im.val[0]);
155 FWD_BOT(x_re.val[1], x_im.val[1], x_re.val[3], x_im.val[3], t_re.val[1], t_im.val[1]);
156 FWD_BOTJ(y_re.val[0], y_im.val[0], y_re.val[2], y_im.val[2], t_re.val[2], t_im.val[2]);
157 FWD_BOTJ(y_re.val[1], y_im.val[1], y_re.val[3], y_im.val[3], t_re.val[3], t_im.val[3]);
158
159 // vstorex4(dre, x_re);
160 // vstorex4(dim, x_im);
161 // vstorex4(dre + 8, y_re);
162 // vstorex4(dim + 8, y_im);
163 // return;
164 }
165 {
166 /*
167 Level 4
168
169 ( 2, 18) * ( 0, 1)
170 ( 3, 19) * ( 0, 1)
171 ( 6, 22) * ( 0, 1)
172 ( 7, 23) * ( 0, 1)
173
174 ( 2, 18) = ( 0, 16) - @
175 ( 3, 19) = ( 1, 17) - @
176 ( 0, 16) = ( 0, 16) + @
177 ( 1, 17) = ( 1, 17) + @
178
179 ( 6, 22) = ( 4, 20) - j@
180 ( 7, 23) = ( 5, 21) - j@
181 ( 4, 20) = ( 4, 20) + j@
182 ( 5, 21) = ( 5, 21) + j@
183
184 ( 10, 26) * ( 2, 3)
185 ( 11, 27) * ( 2, 3)
186 ( 14, 30) * ( 2, 3)
187 ( 15, 31) * ( 2, 3)
188
189 ( 10, 26) = ( 8, 24) - @
190 ( 11, 27) = ( 9, 25) - @
191 ( 8, 24) = ( 8, 24) + @
192 ( 9, 25) = ( 9, 25) + @
193
194 ( 14, 30) = ( 12, 28) - j@
195 ( 15, 31) = ( 13, 29) - j@
196 ( 12, 28) = ( 12, 28) + j@
197 ( 13, 29) = ( 13, 29) + j@
198 */
199 // vloadx4(y_re, dre + 8);
200 // vloadx4(y_im, dim + 8);
201 // vloadx4(x_re, dre);
202 // vloadx4(x_im, dim);
203
204 vloadx2(s_re_im, pom + 4);
205
206 FWD_TOP_LANE(t_re.val[0], t_im.val[0], x_re.val[1], x_im.val[1], s_re_im.val[0]);
207 FWD_TOP_LANE(t_re.val[1], t_im.val[1], x_re.val[3], x_im.val[3], s_re_im.val[0]);
208 FWD_TOP_LANE(t_re.val[2], t_im.val[2], y_re.val[1], y_im.val[1], s_re_im.val[1]);
209 FWD_TOP_LANE(t_re.val[3], t_im.val[3], y_re.val[3], y_im.val[3], s_re_im.val[1]);
210
211 FWD_BOT(x_re.val[0], x_im.val[0], x_re.val[1], x_im.val[1], t_re.val[0], t_im.val[0]);
212 FWD_BOTJ(x_re.val[2], x_im.val[2], x_re.val[3], x_im.val[3], t_re.val[1], t_im.val[1]);
213 FWD_BOT(y_re.val[0], y_im.val[0], y_re.val[1], y_im.val[1], t_re.val[2], t_im.val[2]);
214 FWD_BOTJ(y_re.val[2], y_im.val[2], y_re.val[3], y_im.val[3], t_re.val[3], t_im.val[3]);
215
216 // vstorex4(dre, x_re);
217 // vstorex4(dim, x_im);
218 // vstorex4(dre + 8, y_re);
219 // vstorex4(dim + 8, y_im);
220 // return;
221 }
222 {
223 /*
224 Level 5
225
226 ( 1, 17) * ( 0, 1)
227 ( 5, 21) * ( 2, 3)
228 ------
229 ( 1, 17) = ( 0, 16) - @
230 ( 5, 21) = ( 4, 20) - @
231 ( 0, 16) = ( 0, 16) + @
232 ( 4, 20) = ( 4, 20) + @
233
234 ( 3, 19) * ( 0, 1)
235 ( 7, 23) * ( 2, 3)
236 ------
237 ( 3, 19) = ( 2, 18) - j@
238 ( 7, 23) = ( 6, 22) - j@
239 ( 2, 18) = ( 2, 18) + j@
240 ( 6, 22) = ( 6, 22) + j@
241
242 ( 9, 25) * ( 4, 5)
243 ( 13, 29) * ( 6, 7)
244 ------
245 ( 9, 25) = ( 8, 24) - @
246 ( 13, 29) = ( 12, 28) - @
247 ( 8, 24) = ( 8, 24) + @
248 ( 12, 28) = ( 12, 28) + @
249
250 ( 11, 27) * ( 4, 5)
251 ( 15, 31) * ( 6, 7)
252 ------
253 ( 11, 27) = ( 10, 26) - j@
254 ( 15, 31) = ( 14, 30) - j@
255 ( 10, 26) = ( 10, 26) + j@
256 ( 14, 30) = ( 14, 30) + j@
257
258 before transpose_f64
259 x_re: 0, 1 | 2, 3 | 4, 5 | 6, 7
260 y_re: 8, 9 | 10, 11 | 12, 13 | 14, 15
261 after transpose_f64
262 x_re: 0, 4 | 2, 6 | 1, 5 | 3, 7
263 y_re: 8, 12| 9, 13| 10, 14 | 11, 15
264 after swap
265 x_re: 0, 4 | 1, 5 | 2, 6 | 3, 7
266 y_re: 8, 12| 10, 14 | 9, 13| 11, 15
267 */
268
269 // vloadx4(y_re, dre + 8);
270 // vloadx4(y_im, dim + 8);
271 // vloadx4(x_re, dre);
272 // vloadx4(x_im, dim);
273
274 transpose_f64(x_re, x_re, v_re, 0, 2, 0);
275 transpose_f64(x_re, x_re, v_re, 1, 3, 1);
276 transpose_f64(x_im, x_im, v_im, 0, 2, 0);
277 transpose_f64(x_im, x_im, v_im, 1, 3, 1);
278
279 v_re.val[0] = x_re.val[2];
280 x_re.val[2] = x_re.val[1];
281 x_re.val[1] = v_re.val[0];
282
283 v_im.val[0] = x_im.val[2];
284 x_im.val[2] = x_im.val[1];
285 x_im.val[1] = v_im.val[0];
286
287 transpose_f64(y_re, y_re, v_re, 0, 2, 2);
288 transpose_f64(y_re, y_re, v_re, 1, 3, 3);
289 transpose_f64(y_im, y_im, v_im, 0, 2, 2);
290 transpose_f64(y_im, y_im, v_im, 1, 3, 3);
291
292 v_re.val[0] = y_re.val[2];
293 y_re.val[2] = y_re.val[1];
294 y_re.val[1] = v_re.val[0];
295
296 v_im.val[0] = y_im.val[2];
297 y_im.val[2] = y_im.val[1];
298 y_im.val[1] = v_im.val[0];
299
300 // double pom8[] = {pom[8], pom[12], pom[9], pom[13]};
301 vload2(s_re_im, pom + 8);
302 // vload2(s_re_im, pom8);
303
304 FWD_TOP(t_re.val[0], t_im.val[0], x_re.val[1], x_im.val[1], s_re_im.val[0], s_re_im.val[1]);
305 FWD_TOP(t_re.val[1], t_im.val[1], x_re.val[3], x_im.val[3], s_re_im.val[0], s_re_im.val[1]);
306
307 // double pom12[] = {pom[10], pom[14], pom[11], pom[15]};
308 vload2(s_re_im, pom + 12);
309 // vload2(s_re_im, pom12);
310
311 FWD_TOP(t_re.val[2], t_im.val[2], y_re.val[1], y_im.val[1], s_re_im.val[0], s_re_im.val[1]);
312 FWD_TOP(t_re.val[3], t_im.val[3], y_re.val[3], y_im.val[3], s_re_im.val[0], s_re_im.val[1]);
313
314 FWD_BOT(x_re.val[0], x_im.val[0], x_re.val[1], x_im.val[1], t_re.val[0], t_im.val[0]);
315 FWD_BOTJ(x_re.val[2], x_im.val[2], x_re.val[3], x_im.val[3], t_re.val[1], t_im.val[1]);
316
317 vstore4(dre, x_re);
318 vstore4(dim, x_im);
319
320 FWD_BOT(y_re.val[0], y_im.val[0], y_re.val[1], y_im.val[1], t_re.val[2], t_im.val[2]);
321 FWD_BOTJ(y_re.val[2], y_im.val[2], y_re.val[3], y_im.val[3], t_re.val[3], t_im.val[3]);
322
323 vstore4(dre + 8, y_re);
324 vstore4(dim + 8, y_im);
325 }
326}
327
328void reim_twiddle_fft_neon(uint64_t h, double* re, double* im, double om[2]) {
329 // Total SIMD register: 28 = 24 + 4
330 if (h < 8) abort(); // bug
331 float64x2_t s_re_im; // 2
332 float64x2x4_t x_re, x_im, y_re, y_im, v_re, v_im; // 32
333 vload(s_re_im, om);
334 for (uint64_t blk = 0; blk < h; blk += 8) {
335 double* dre = re + blk;
336 double* dim = im + blk;
337 vloadx4(y_re, dre + h);
338 vloadx4(y_im, dim + h);
339 FWD_TOP_LANEx4(v_re, v_im, y_re, y_im, s_re_im);
340 vloadx4(x_re, dre);
341 vloadx4(x_im, dim);
342 FWD_BOTx4(x_re, x_im, y_re, y_im, v_re, v_im);
343 vstorex4(dre, x_re);
344 vstorex4(dim, x_im);
345 vstorex4(dre + h, y_re);
346 vstorex4(dim + h, y_im);
347 }
348}
349
350void reim_ctwiddle(double* ra, double* ia, double* rb, double* ib, double omre, double omim);
351// i (omre + i omim) = -omim + i omre
352void reim_citwiddle(double* ra, double* ia, double* rb, double* ib, double omre, double omim);
353
354void reim_bitwiddle_fft_neon(uint64_t h, double* re, double* im, double om[4]) {
355 // Total SIMD register: 28 = 24 + 4
356 if (h < 4) abort(); // bug
357 double* r0 = re;
358 double* r1 = re + h;
359 double* r2 = re + 2 * h;
360 double* r3 = re + 3 * h;
361 double* i0 = im;
362 double* i1 = im + h;
363 double* i2 = im + 2 * h;
364 double* i3 = im + 3 * h;
365 float64x2x2_t s_re_im; // 2
366 float64x2x4_t v_re, v_im; // 2
367 float64x2x2_t x0_re, x0_im, x1_re, x1_im; // 32
368 float64x2x2_t x2_re, x2_im, x3_re, x3_im; // 32
369 vloadx2(s_re_im, om);
370 for (uint64_t blk = 0; blk < h; blk += 4) {
371 {
372 vloadx2(x2_re, r2 + blk);
373 vloadx2(x3_re, r3 + blk);
374 vloadx2(x2_im, i2 + blk);
375 vloadx2(x3_im, i3 + blk);
376 FWD_TOP_LANE(v_re.val[0], v_im.val[0], x2_re.val[0], x2_im.val[0], s_re_im.val[0]);
377 FWD_TOP_LANE(v_re.val[1], v_im.val[1], x2_re.val[1], x2_im.val[1], s_re_im.val[0]);
378 FWD_TOP_LANE(v_re.val[2], v_im.val[2], x3_re.val[0], x3_im.val[0], s_re_im.val[0]);
379 FWD_TOP_LANE(v_re.val[3], v_im.val[3], x3_re.val[1], x3_im.val[1], s_re_im.val[0]);
380 vloadx2(x0_re, r0 + blk);
381 vloadx2(x1_re, r1 + blk);
382 vloadx2(x0_im, i0 + blk);
383 vloadx2(x1_im, i1 + blk);
384 FWD_BOT(x0_re.val[0], x0_im.val[0], x2_re.val[0], x2_im.val[0], v_re.val[0], v_im.val[0]);
385 FWD_BOT(x0_re.val[1], x0_im.val[1], x2_re.val[1], x2_im.val[1], v_re.val[1], v_im.val[1]);
386 FWD_BOT(x1_re.val[0], x1_im.val[0], x3_re.val[0], x3_im.val[0], v_re.val[2], v_im.val[2]);
387 FWD_BOT(x1_re.val[1], x1_im.val[1], x3_re.val[1], x3_im.val[1], v_re.val[3], v_im.val[3]);
388 }
389 {
390 FWD_TOP_LANE(v_re.val[0], v_im.val[0], x1_re.val[0], x1_im.val[0], s_re_im.val[1]);
391 FWD_TOP_LANE(v_re.val[1], v_im.val[1], x1_re.val[1], x1_im.val[1], s_re_im.val[1]);
392 FWD_TOP_LANE(v_re.val[2], v_im.val[2], x3_re.val[0], x3_im.val[0], s_re_im.val[1]);
393 FWD_TOP_LANE(v_re.val[3], v_im.val[3], x3_re.val[1], x3_im.val[1], s_re_im.val[1]);
394 FWD_BOT(x0_re.val[0], x0_im.val[0], x1_re.val[0], x1_im.val[0], v_re.val[0], v_im.val[0]);
395 FWD_BOT(x0_re.val[1], x0_im.val[1], x1_re.val[1], x1_im.val[1], v_re.val[1], v_im.val[1]);
396 FWD_BOTJ(x2_re.val[0], x2_im.val[0], x3_re.val[0], x3_im.val[0], v_re.val[2], v_im.val[2]);
397 FWD_BOTJ(x2_re.val[1], x2_im.val[1], x3_re.val[1], x3_im.val[1], v_re.val[3], v_im.val[3]);
398 vstorex2(r0 + blk, x0_re);
399 vstorex2(r1 + blk, x1_re);
400 vstorex2(r2 + blk, x2_re);
401 vstorex2(r3 + blk, x3_re);
402 vstorex2(i0 + blk, x0_im);
403 vstorex2(i1 + blk, x1_im);
404 vstorex2(i2 + blk, x2_im);
405 vstorex2(i3 + blk, x3_im);
406 }
407 }
408}
409
410#if 0
411static void ZfN(iFFT_log2)(fpr *f)
412{
413 /*
414 y_re: 1 = (2 - 3) * 5 + (0 - 1) * 4
415 y_im: 3 = (2 - 3) * 4 - (0 - 1) * 5
416 x_re: 0 = 0 + 1
417 x_im: 2 = 2 + 3
418
419 Turn out this vectorize code is too short to be executed,
420 the scalar version is consistently faster
421
422 float64x2x2_t tmp;
423 float64x2_t v, s, t;
424
425 // 0: 0, 2
426 // 1: 1, 3
427
428 vload2(tmp, &f[0]);
429 vload(s, &fpr_gm_tab[4]);
430
431 vfsub(v, tmp.val[0], tmp.val[1]);
432 vfadd(tmp.val[0], tmp.val[0], tmp.val[1]);
433
434 // y_im: 3 = (2 - 3) * 4 - (0 - 1) * 5
435 // y_re: 1 = (2 - 3) * 5 + (0 - 1) * 4
436 vswap(t, v);
437
438 vfmul_lane(tmp.val[1], s, t, 0);
439 vfcmla_90(tmp.val[1], t, s);
440
441 vfmuln(tmp.val[0], tmp.val[0], 0.5);
442 vfmuln(tmp.val[1], tmp.val[1], 0.5);
443
444 vswap(tmp.val[1], tmp.val[1]);
445
446 vstore2(&f[0], tmp);
447 */
448
449 fpr x_re, x_im, y_re, y_im, s;
450 x_re = f[0];
451 y_re = f[1];
452 x_im = f[2];
453 y_im = f[3];
454 s = fpr_tab_log2[0] * 0.5;
455
456 f[0] = (x_re + y_re) * 0.5;
457 f[2] = (x_im + y_im) * 0.5;
458
459 x_re = (x_re - y_re) * s;
460 x_im = (x_im - y_im) * s;
461
462 f[1] = x_im + x_re;
463 f[3] = x_im - x_re;
464}
465
466static void ZfN(iFFT_log3)(fpr *f)
467{
468 /*
469 * Total instructions: 27
470 y_re: 1 = (4 - 5) * 9 + (0 - 1) * 8
471 y_re: 3 = (6 - 7) * 11 + (2 - 3) * 10
472 y_im: 5 = (4 - 5) * 8 - (0 - 1) * 9
473 y_im: 7 = (6 - 7) * 10 - (2 - 3) * 11
474 x_re: 0 = 0 + 1
475 x_re: 2 = 2 + 3
476 x_im: 4 = 4 + 5
477 x_im: 6 = 6 + 7
478 */
479 // 0: 0, 2 - 0: 0, 4
480 // 1: 1, 3 - 1: 1, 5
481 // 2: 4, 6 - 2: 2, 6
482 // 3: 5, 7 - 3: 3, 7
483 float64x2x4_t tmp;
484 float64x2x2_t x_re_im, y_re_im, v, s_re_im;
485
486 vload2(x_re_im, &f[0]);
487 vload2(y_re_im, &f[4]);
488
489 vfsub(v.val[0], x_re_im.val[0], x_re_im.val[1]);
490 vfsub(v.val[1], y_re_im.val[0], y_re_im.val[1]);
491 vfadd(x_re_im.val[0], x_re_im.val[0], x_re_im.val[1]);
492 vfadd(x_re_im.val[1], y_re_im.val[0], y_re_im.val[1]);
493
494 // 0: 8, 10
495 // 1: 9, 11
496 vload2(s_re_im, &fpr_tab_log3[0]);
497
498 vfmul(y_re_im.val[0], v.val[1], s_re_im.val[1]);
499 vfmla(y_re_im.val[0], y_re_im.val[0], v.val[0], s_re_im.val[0]);
500 vfmul(y_re_im.val[1], v.val[1], s_re_im.val[0]);
501 vfmls(y_re_im.val[1], y_re_im.val[1], v.val[0], s_re_im.val[1]);
502
503 // x: 0,2 | 4,6
504 // y: 1,3 | 5,7
505 tmp.val[0] = vtrn1q_f64(x_re_im.val[0], y_re_im.val[0]);
506 tmp.val[1] = vtrn2q_f64(x_re_im.val[0], y_re_im.val[0]);
507 tmp.val[2] = vtrn1q_f64(x_re_im.val[1], y_re_im.val[1]);
508 tmp.val[3] = vtrn2q_f64(x_re_im.val[1], y_re_im.val[1]);
509 // tmp: 0,1 | 2,3 | 4,5 | 6,7
510 /*
511 y_re: 2 = (4 - 6) * 4 + (0 - 2) * 4
512 y_re: 3 = (5 - 7) * 4 + (1 - 3) * 4
513 y_im: 6 = (4 - 6) * 4 - (0 - 2) * 4
514 y_im: 7 = (5 - 7) * 4 - (1 - 3) * 4
515 x_re: 0 = 0 + 2
516 x_re: 1 = 1 + 3
517 x_im: 4 = 4 + 6
518 x_im: 5 = 5 + 7
519 */
520 s_re_im.val[0] = vld1q_dup_f64(&fpr_tab_log2[0]);
521
522 vfadd(x_re_im.val[0], tmp.val[0], tmp.val[1]);
523 vfadd(x_re_im.val[1], tmp.val[2], tmp.val[3]);
524 vfsub(v.val[0], tmp.val[0], tmp.val[1]);
525 vfsub(v.val[1], tmp.val[2], tmp.val[3]);
526
527 vfmuln(tmp.val[0], x_re_im.val[0], 0.25);
528 vfmuln(tmp.val[2], x_re_im.val[1], 0.25);
529
530 vfmuln(s_re_im.val[0], s_re_im.val[0], 0.25);
531
532 vfmul(y_re_im.val[0], v.val[0], s_re_im.val[0]);
533 vfmul(y_re_im.val[1], v.val[1], s_re_im.val[0]);
534
535 vfadd(tmp.val[1], y_re_im.val[1], y_re_im.val[0]);
536 vfsub(tmp.val[3], y_re_im.val[1], y_re_im.val[0]);
537
538 vstorex4(&f[0], tmp);
539}
540
541static void ZfN(iFFT_log4)(fpr *f)
542{
543 /*
544 * ( 0, 8) - ( 1, 9)
545 * ( 4, 12) - ( 5, 13)
546 *
547 * ( 0, 8) + ( 1, 9)
548 * ( 4, 12) + ( 5, 13)
549 *
550 * ( 3, 11) - ( 2, 10)
551 * ( 7, 15) - ( 6, 14)
552 *
553 * ( 2, 10) + ( 3, 11)
554 * ( 6, 14) + ( 7, 15)
555 *
556 * ( 1, 9) = @ * ( 0, 1)
557 * ( 5, 13) = @ * ( 2, 3)
558 *
559 * ( 3, 11) = j@ * ( 0, 1)
560 * ( 7, 15) = j@ * ( 2, 3)
561 */
562
563 float64x2x4_t re, im, t;
564 float64x2x2_t t_re, t_im, s_re_im;
565
566 vload4(re, &f[0]);
567 vload4(im, &f[8]);
568
569 INV_TOPJ (t_re.val[0], t_im.val[0], re.val[0], im.val[0], re.val[1], im.val[1]);
570 INV_TOPJm(t_re.val[1], t_im.val[1], re.val[2], im.val[2], re.val[3], im.val[3]);
571
572 vload2(s_re_im, &fpr_tab_log4[0]);
573
574 INV_BOTJ (re.val[1], im.val[1], t_re.val[0], t_im.val[0], s_re_im.val[0], s_re_im.val[1]);
575 INV_BOTJm(re.val[3], im.val[3], t_re.val[1], t_im.val[1], s_re_im.val[0], s_re_im.val[1]);
576
577 /*
578 * ( 0, 8) - ( 2, 10)
579 * ( 1, 9) - ( 3, 11)
580 *
581 * ( 0, 8) + ( 2, 10)
582 * ( 1, 9) + ( 3, 11)
583 *
584 * ( 2, 10) = @ * ( 0, 1)
585 * ( 3, 11) = @ * ( 0, 1)
586 *
587 * ( 6, 14) - ( 4, 12)
588 * ( 7, 15) - ( 5, 13)
589 *
590 * ( 4, 12) + ( 6, 14)
591 * ( 5, 13) + ( 7, 15)
592 *
593 * ( 6, 14) = j@ * ( 0, 1)
594 * ( 7, 15) = j@ * ( 0, 1)
595 */
596
597 // re: 0, 4 | 1, 5 | 2, 6 | 3, 7
598 // im: 8, 12| 9, 13|10, 14|11, 15
599
600 transpose_f64(re, re, t, 0, 1, 0);
601 transpose_f64(re, re, t, 2, 3, 1);
602 transpose_f64(im, im, t, 0, 1, 2);
603 transpose_f64(im, im, t, 2, 3, 3);
604
605 // re: 0, 1 | 4, 5 | 2, 3 | 6, 7
606 // im: 8, 9 | 12, 13|10, 11| 14, 15
607 t.val[0] = re.val[1];
608 re.val[1] = re.val[2];
609 re.val[2] = t.val[0];
610
611 t.val[1] = im.val[1];
612 im.val[1] = im.val[2];
613 im.val[2] = t.val[1];
614
615 // re: 0, 1 | 2, 3| 4, 5 | 6, 7
616 // im: 8, 9 | 10, 11| 12, 13| 14, 15
617
618 INV_TOPJ (t_re.val[0], t_im.val[0], re.val[0], im.val[0], re.val[1], im.val[1]);
619 INV_TOPJm(t_re.val[1], t_im.val[1], re.val[2], im.val[2], re.val[3], im.val[3]);
620
621 vload(s_re_im.val[0], &fpr_tab_log3[0]);
622
623 INV_BOTJ_LANE (re.val[1], im.val[1], t_re.val[0], t_im.val[0], s_re_im.val[0]);
624 INV_BOTJm_LANE(re.val[3], im.val[3], t_re.val[1], t_im.val[1], s_re_im.val[0]);
625
626 /*
627 * ( 0, 8) - ( 4, 12)
628 * ( 1, 9) - ( 5, 13)
629 * ( 0, 8) + ( 4, 12)
630 * ( 1, 9) + ( 5, 13)
631 *
632 * ( 2, 10) - ( 6, 14)
633 * ( 3, 11) - ( 7, 15)
634 * ( 2, 10) + ( 6, 14)
635 * ( 3, 11) + ( 7, 15)
636 *
637 * ( 4, 12) = @ * ( 0, 1)
638 * ( 5, 13) = @ * ( 0, 1)
639 *
640 * ( 6, 14) = @ * ( 0, 1)
641 * ( 7, 15) = @ * ( 0, 1)
642 */
643
644 INV_TOPJ(t_re.val[0], t_im.val[0], re.val[0], im.val[0], re.val[2], im.val[2]);
645 INV_TOPJ(t_re.val[1], t_im.val[1], re.val[1], im.val[1], re.val[3], im.val[3]);
646
647 vfmuln(re.val[0], re.val[0], 0.12500000000);
648 vfmuln(re.val[1], re.val[1], 0.12500000000);
649 vfmuln(im.val[0], im.val[0], 0.12500000000);
650 vfmuln(im.val[1], im.val[1], 0.12500000000);
651
652 s_re_im.val[0] = vld1q_dup_f64(&fpr_tab_log2[0]);
653
654 vfmuln(s_re_im.val[0], s_re_im.val[0], 0.12500000000);
655
656 vfmul(t_re.val[0], t_re.val[0], s_re_im.val[0]);
657 vfmul(t_re.val[1], t_re.val[1], s_re_im.val[0]);
658 vfmul(t_im.val[0], t_im.val[0], s_re_im.val[0]);
659 vfmul(t_im.val[1], t_im.val[1], s_re_im.val[0]);
660
661 vfsub(im.val[2], t_im.val[0], t_re.val[0]);
662 vfsub(im.val[3], t_im.val[1], t_re.val[1]);
663 vfadd(re.val[2], t_im.val[0], t_re.val[0]);
664 vfadd(re.val[3], t_im.val[1], t_re.val[1]);
665
666 vstorex4(&f[0], re);
667 vstorex4(&f[8], im);
668}
669
670static
671 void ZfN(iFFT_log5)(fpr *f, const unsigned logn, const unsigned last)
672{
673 // Total SIMD register: 26 = 24 + 2
674 float64x2x4_t x_re, x_im, y_re, y_im, t_re, t_im; // 24
675 float64x2x2_t s_re_im; // 2
676 const unsigned n = 1 << logn;
677 const unsigned hn = n >> 1;
678
679 int level = logn;
680 const fpr *fpr_tab5 = fpr_table[level--],
681 *fpr_tab4 = fpr_table[level--],
682 *fpr_tab3 = fpr_table[level--],
683 *fpr_tab2 = fpr_table[level];
684 int k2 = 0, k3 = 0, k4 = 0, k5 = 0;
685
686 for (unsigned j = 0; j < hn; j += 16)
687 {
688 /*
689 * ( 0, 16) - ( 1, 17)
690 * ( 4, 20) - ( 5, 21)
691 * ( 0, 16) + ( 1, 17)
692 * ( 4, 20) + ( 5, 21)
693 * ( 1, 17) = @ * ( 0, 1)
694 * ( 5, 21) = @ * ( 2, 3)
695 *
696 * ( 2, 18) - ( 3, 19)
697 * ( 6, 22) - ( 7, 23)
698 * ( 2, 18) + ( 3, 19)
699 * ( 6, 22) + ( 7, 23)
700 * ( 3, 19) = j@ * ( 0, 1)
701 * ( 7, 23) = j@ * ( 2, 3)
702 *
703 * ( 8, 24) - ( 9, 25)
704 * ( 12, 28) - ( 13, 29)
705 * ( 8, 24) + ( 9, 25)
706 * ( 12, 28) + ( 13, 29)
707 * ( 9, 25) = @ * ( 4, 5)
708 * ( 13, 29) = @ * ( 6, 7)
709 *
710 * ( 10, 26) - ( 11, 27)
711 * ( 14, 30) - ( 15, 31)
712 * ( 10, 26) + ( 11, 27)
713 * ( 14, 30) + ( 15, 31)
714 * ( 11, 27) = j@ * ( 4, 5)
715 * ( 15, 31) = j@ * ( 6, 7)
716 */
717
718 vload4(x_re, &f[j]);
719 vload4(x_im, &f[j + hn]);
720
721 INV_TOPJ(t_re.val[0], t_im.val[0], x_re.val[0], x_im.val[0], x_re.val[1], x_im.val[1]);
722 INV_TOPJm(t_re.val[2], t_im.val[2], x_re.val[2], x_im.val[2], x_re.val[3], x_im.val[3]);
723
724 vload4(y_re, &f[j + 8]);
725 vload4(y_im, &f[j + 8 + hn])
726
727 INV_TOPJ(t_re.val[1], t_im.val[1], y_re.val[0], y_im.val[0], y_re.val[1], y_im.val[1]);
728 INV_TOPJm(t_re.val[3], t_im.val[3], y_re.val[2], y_im.val[2], y_re.val[3], y_im.val[3]);
729
730 vload2(s_re_im, &fpr_tab5[k5]);
731 k5 += 4;
732
733 INV_BOTJ (x_re.val[1], x_im.val[1], t_re.val[0], t_im.val[0], s_re_im.val[0], s_re_im.val[1]);
734 INV_BOTJm(x_re.val[3], x_im.val[3], t_re.val[2], t_im.val[2], s_re_im.val[0], s_re_im.val[1]);
735
736 vload2(s_re_im, &fpr_tab5[k5]);
737 k5 += 4;
738
739 INV_BOTJ (y_re.val[1], y_im.val[1], t_re.val[1], t_im.val[1], s_re_im.val[0], s_re_im.val[1]);
740 INV_BOTJm(y_re.val[3], y_im.val[3], t_re.val[3], t_im.val[3], s_re_im.val[0], s_re_im.val[1]);
741
742
743 // x_re: 0, 4 | 1, 5 | 2, 6 | 3, 7
744 // y_re: 8, 12| 9, 13|10, 14|11, 15
745
746 transpose_f64(x_re, x_re, t_re, 0, 1, 0);
747 transpose_f64(x_re, x_re, t_re, 2, 3, 1);
748 transpose_f64(y_re, y_re, t_re, 0, 1, 2);
749 transpose_f64(y_re, y_re, t_re, 2, 3, 3);
750
751 transpose_f64(x_im, x_im, t_im, 0, 1, 0);
752 transpose_f64(x_im, x_im, t_im, 2, 3, 1);
753 transpose_f64(y_im, y_im, t_im, 0, 1, 2);
754 transpose_f64(y_im, y_im, t_im, 2, 3, 3);
755
756 // x_re: 0, 1 | 4, 5 | 2, 3 | 6, 7
757 // y_re: 8, 9 | 12,13|10,11 |14, 15
758
759 t_re.val[0] = x_re.val[1];
760 x_re.val[1] = x_re.val[2];
761 x_re.val[2] = t_re.val[0];
762
763 t_re.val[1] = y_re.val[1];
764 y_re.val[1] = y_re.val[2];
765 y_re.val[2] = t_re.val[1];
766
767
768 t_im.val[0] = x_im.val[1];
769 x_im.val[1] = x_im.val[2];
770 x_im.val[2] = t_im.val[0];
771
772 t_im.val[1] = y_im.val[1];
773 y_im.val[1] = y_im.val[2];
774 y_im.val[2] = t_im.val[1];
775 // x_re: 0, 1 | 2, 3| 4, 5 | 6, 7
776 // y_re: 8, 9 | 10, 11| 12, 13| 14, 15
777
778 /*
779 * ( 0, 16) - ( 2, 18)
780 * ( 1, 17) - ( 3, 19)
781 * ( 0, 16) + ( 2, 18)
782 * ( 1, 17) + ( 3, 19)
783 * ( 2, 18) = @ * ( 0, 1)
784 * ( 3, 19) = @ * ( 0, 1)
785 *
786 * ( 4, 20) - ( 6, 22)
787 * ( 5, 21) - ( 7, 23)
788 * ( 4, 20) + ( 6, 22)
789 * ( 5, 21) + ( 7, 23)
790 * ( 6, 22) = j@ * ( 0, 1)
791 * ( 7, 23) = j@ * ( 0, 1)
792 *
793 * ( 8, 24) - ( 10, 26)
794 * ( 9, 25) - ( 11, 27)
795 * ( 8, 24) + ( 10, 26)
796 * ( 9, 25) + ( 11, 27)
797 * ( 10, 26) = @ * ( 2, 3)
798 * ( 11, 27) = @ * ( 2, 3)
799 *
800 * ( 12, 28) - ( 14, 30)
801 * ( 13, 29) - ( 15, 31)
802 * ( 12, 28) + ( 14, 30)
803 * ( 13, 29) + ( 15, 31)
804 * ( 14, 30) = j@ * ( 2, 3)
805 * ( 15, 31) = j@ * ( 2, 3)
806 */
807
808 INV_TOPJ (t_re.val[0], t_im.val[0], x_re.val[0], x_im.val[0], x_re.val[1], x_im.val[1]);
809 INV_TOPJm(t_re.val[1], t_im.val[1], x_re.val[2], x_im.val[2], x_re.val[3], x_im.val[3]);
810
811 INV_TOPJ (t_re.val[2], t_im.val[2], y_re.val[0], y_im.val[0], y_re.val[1], y_im.val[1]);
812 INV_TOPJm(t_re.val[3], t_im.val[3], y_re.val[2], y_im.val[2], y_re.val[3], y_im.val[3]);
813
814 vloadx2(s_re_im, &fpr_tab4[k4]);
815 k4 += 4;
816
817 INV_BOTJ_LANE (x_re.val[1], x_im.val[1], t_re.val[0], t_im.val[0], s_re_im.val[0]);
818 INV_BOTJm_LANE(x_re.val[3], x_im.val[3], t_re.val[1], t_im.val[1], s_re_im.val[0]);
819
820 INV_BOTJ_LANE (y_re.val[1], y_im.val[1], t_re.val[2], t_im.val[2], s_re_im.val[1]);
821 INV_BOTJm_LANE(y_re.val[3], y_im.val[3], t_re.val[3], t_im.val[3], s_re_im.val[1]);
822
823 /*
824 * ( 0, 16) - ( 4, 20)
825 * ( 1, 17) - ( 5, 21)
826 * ( 0, 16) + ( 4, 20)
827 * ( 1, 17) + ( 5, 21)
828 * ( 4, 20) = @ * ( 0, 1)
829 * ( 5, 21) = @ * ( 0, 1)
830 *
831 * ( 2, 18) - ( 6, 22)
832 * ( 3, 19) - ( 7, 23)
833 * ( 2, 18) + ( 6, 22)
834 * ( 3, 19) + ( 7, 23)
835 * ( 6, 22) = @ * ( 0, 1)
836 * ( 7, 23) = @ * ( 0, 1)
837 *
838 * ( 8, 24) - ( 12, 28)
839 * ( 9, 25) - ( 13, 29)
840 * ( 8, 24) + ( 12, 28)
841 * ( 9, 25) + ( 13, 29)
842 * ( 12, 28) = j@ * ( 0, 1)
843 * ( 13, 29) = j@ * ( 0, 1)
844 *
845 * ( 10, 26) - ( 14, 30)
846 * ( 11, 27) - ( 15, 31)
847 * ( 10, 26) + ( 14, 30)
848 * ( 11, 27) + ( 15, 31)
849 * ( 14, 30) = j@ * ( 0, 1)
850 * ( 15, 31) = j@ * ( 0, 1)
851 */
852
853 INV_TOPJ (t_re.val[0], t_im.val[0], x_re.val[0], x_im.val[0], x_re.val[2], x_im.val[2]);
854 INV_TOPJ (t_re.val[1], t_im.val[1], x_re.val[1], x_im.val[1], x_re.val[3], x_im.val[3]);
855
856 INV_TOPJm(t_re.val[2], t_im.val[2], y_re.val[0], y_im.val[0], y_re.val[2], y_im.val[2]);
857 INV_TOPJm(t_re.val[3], t_im.val[3], y_re.val[1], y_im.val[1], y_re.val[3], y_im.val[3]);
858
859 vload(s_re_im.val[0], &fpr_tab3[k3]);
860 k3 += 2;
861
862 INV_BOTJ_LANE(x_re.val[2], x_im.val[2], t_re.val[0], t_im.val[0], s_re_im.val[0]);
863 INV_BOTJ_LANE(x_re.val[3], x_im.val[3], t_re.val[1], t_im.val[1], s_re_im.val[0]);
864
865 INV_BOTJm_LANE(y_re.val[2], y_im.val[2], t_re.val[2], t_im.val[2], s_re_im.val[0]);
866 INV_BOTJm_LANE(y_re.val[3], y_im.val[3], t_re.val[3], t_im.val[3], s_re_im.val[0]);
867
868 /*
869 * ( 0, 16) - ( 8, 24)
870 * ( 1, 17) - ( 9, 25)
871 * ( 0, 16) + ( 8, 24)
872 * ( 1, 17) + ( 9, 25)
873 * ( 8, 24) = @ * ( 0, 1)
874 * ( 9, 25) = @ * ( 0, 1)
875 *
876 * ( 2, 18) - ( 10, 26)
877 * ( 3, 19) - ( 11, 27)
878 * ( 2, 18) + ( 10, 26)
879 * ( 3, 19) + ( 11, 27)
880 * ( 10, 26) = @ * ( 0, 1)
881 * ( 11, 27) = @ * ( 0, 1)
882 *
883 * ( 4, 20) - ( 12, 28)
884 * ( 5, 21) - ( 13, 29)
885 * ( 4, 20) + ( 12, 28)
886 * ( 5, 21) + ( 13, 29)
887 * ( 12, 28) = @ * ( 0, 1)
888 * ( 13, 29) = @ * ( 0, 1)
889 *
890 * ( 6, 22) - ( 14, 30)
891 * ( 7, 23) - ( 15, 31)
892 * ( 6, 22) + ( 14, 30)
893 * ( 7, 23) + ( 15, 31)
894 * ( 14, 30) = @ * ( 0, 1)
895 * ( 15, 31) = @ * ( 0, 1)
896 */
897
898
899 if ( (j >> 4) & 1)
900 {
901 INV_TOPJmx4(t_re, t_im, x_re, x_im, y_re, y_im);
902 }
903 else
904 {
905 INV_TOPJx4(t_re, t_im, x_re, x_im, y_re, y_im);
906 }
907
908 vload(s_re_im.val[0], &fpr_tab2[k2]);
909 k2 += 2 * ((j & 31) == 16);
910
911 if (last)
912 {
913 vfmuln(s_re_im.val[0], s_re_im.val[0], fpr_p2_tab[logn]);
914 vfmulnx4(x_re, x_re, fpr_p2_tab[logn]);
915 vfmulnx4(x_im, x_im, fpr_p2_tab[logn]);
916 }
917 vstorex4(&f[j], x_re);
918 vstorex4(&f[j + hn], x_im);
919
920 if (logn == 5)
921 {
922 // Special case in fpr_tab_log2 where re == im
923 vfmulx4_i(t_re, t_re, s_re_im.val[0]);
924 vfmulx4_i(t_im, t_im, s_re_im.val[0]);
925
926 vfaddx4(y_re, t_im, t_re);
927 vfsubx4(y_im, t_im, t_re);
928 }
929 else
930 {
931 if ((j >> 4) & 1)
932 {
933 INV_BOTJm_LANEx4(y_re, y_im, t_re, t_im, s_re_im.val[0]);
934 }
935 else
936 {
937 INV_BOTJ_LANEx4(y_re, y_im, t_re, t_im, s_re_im.val[0]);
938 }
939 }
940
941 vstorex4(&f[j + 8], y_re);
942 vstorex4(&f[j + 8 + hn], y_im);
943 }
944}
945
946static
947 void ZfN(iFFT_logn1)(fpr *f, const unsigned logn, const unsigned last)
948{
949 // Total SIMD register 26 = 24 + 2
950 float64x2x4_t a_re, a_im, b_re, b_im, t_re, t_im; // 24
951 float64x2_t s_re_im; // 2
952
953 const unsigned n = 1 << logn;
954 const unsigned hn = n >> 1;
955 const unsigned ht = n >> 2;
956
957 for (unsigned j = 0; j < ht; j+= 8)
958 {
959 vloadx4(a_re, &f[j]);
960 vloadx4(a_im, &f[j + hn]);
961 vloadx4(b_re, &f[j + ht]);
962 vloadx4(b_im, &f[j + ht + hn]);
963
964 INV_TOPJx4(t_re, t_im, a_re, a_im, b_re, b_im);
965
966 s_re_im = vld1q_dup_f64(&fpr_tab_log2[0]);
967
968 if (last)
969 {
970 vfmuln(s_re_im, s_re_im, fpr_p2_tab[logn]);
971 vfmulnx4(a_re, a_re, fpr_p2_tab[logn]);
972 vfmulnx4(a_im, a_im, fpr_p2_tab[logn]);
973 }
974
975 vstorex4(&f[j], a_re);
976 vstorex4(&f[j + hn], a_im);
977
978 vfmulx4_i(t_re, t_re, s_re_im);
979 vfmulx4_i(t_im, t_im, s_re_im);
980
981 vfaddx4(b_re, t_im, t_re);
982 vfsubx4(b_im, t_im, t_re);
983
984 vstorex4(&f[j + ht], b_re);
985 vstorex4(&f[j + ht + hn], b_im);
986 }
987}
988
989// static
990// void ZfN(iFFT_logn2)(fpr *f, const unsigned logn, const unsigned level, unsigned last)
991// {
992// const unsigned int falcon_n = 1 << logn;
993// const unsigned int hn = falcon_n >> 1;
994
995// // Total SIMD register: 26 = 16 + 8 + 2
996// float64x2x4_t t_re, t_im; // 8
997// float64x2x2_t x1_re, x2_re, x1_im, x2_im,
998// y1_re, y2_re, y1_im, y2_im; // 16
999// float64x2_t s1_re_im, s2_re_im; // 2
1000
1001// const fpr *fpr_inv_tab1 = NULL, *fpr_inv_tab2 = NULL;
1002// unsigned l, len, start, j, k1, k2;
1003// unsigned bar = logn - 4 - 2;
1004// unsigned Jm;
1005
1006// for (l = 4; l < logn - level - 1; l += 2)
1007// {
1008// len = 1 << l;
1009// last -= 1;
1010// fpr_inv_tab1 = fpr_table[bar--];
1011// fpr_inv_tab2 = fpr_table[bar--];
1012// k1 = 0; k2 = 0;
1013
1014// for (start = 0; start < hn; start += 1 << (l + 2))
1015// {
1016// vload(s1_re_im, &fpr_inv_tab1[k1]);
1017// vload(s2_re_im, &fpr_inv_tab2[k2]);
1018// k1 += 2;
1019// k2 += 2 * ((start & 127) == 64);
1020// if (!last)
1021// {
1022// vfmuln(s2_re_im, s2_re_im, fpr_p2_tab[logn]);
1023// }
1024// Jm = (start >> (l+ 2)) & 1;
1025// for (j = start; j < start + len; j += 4)
1026// {
1027// /*
1028// Level 6
1029// * ( 0, 64) - ( 16, 80)
1030// * ( 1, 65) - ( 17, 81)
1031// * ( 0, 64) + ( 16, 80)
1032// * ( 1, 65) + ( 17, 81)
1033// * ( 16, 80) = @ * ( 0, 1)
1034// * ( 17, 81) = @ * ( 0, 1)
1035// *
1036// * ( 2, 66) - ( 18, 82)
1037// * ( 3, 67) - ( 19, 83)
1038// * ( 2, 66) + ( 18, 82)
1039// * ( 3, 67) + ( 19, 83)
1040// * ( 18, 82) = @ * ( 0, 1)
1041// * ( 19, 83) = @ * ( 0, 1)
1042// *
1043// * ( 32, 96) - ( 48, 112)
1044// * ( 33, 97) - ( 49, 113)
1045// * ( 32, 96) + ( 48, 112)
1046// * ( 33, 97) + ( 49, 113)
1047// * ( 48, 112) = j@ * ( 0, 1)
1048// * ( 49, 113) = j@ * ( 0, 1)
1049// *
1050// * ( 34, 98) - ( 50, 114)
1051// * ( 35, 99) - ( 51, 115)
1052// * ( 34, 98) + ( 50, 114)
1053// * ( 35, 99) + ( 51, 115)
1054// * ( 50, 114) = j@ * ( 0, 1)
1055// * ( 51, 115) = j@ * ( 0, 1)
1056// */
1057// // x1: 0 -> 4 | 64 -> 67
1058// // y1: 16 -> 19 | 80 -> 81
1059// // x2: 32 -> 35 | 96 -> 99
1060// // y2: 48 -> 51 | 112 -> 115
1061// vloadx2(x1_re, &f[j]);
1062// vloadx2(x1_im, &f[j + hn]);
1063// vloadx2(y1_re, &f[j + len]);
1064// vloadx2(y1_im, &f[j + len + hn]);
1065
1066// INV_TOPJ (t_re.val[0], t_im.val[0], x1_re.val[0], x1_im.val[0], y1_re.val[0], y1_im.val[0]);
1067// INV_TOPJ (t_re.val[1], t_im.val[1], x1_re.val[1], x1_im.val[1], y1_re.val[1], y1_im.val[1]);
1068
1069// vloadx2(x2_re, &f[j + 2*len]);
1070// vloadx2(x2_im, &f[j + 2*len + hn]);
1071// vloadx2(y2_re, &f[j + 3*len]);
1072// vloadx2(y2_im, &f[j + 3*len + hn]);
1073
1074// INV_TOPJm(t_re.val[2], t_im.val[2], x2_re.val[0], x2_im.val[0], y2_re.val[0], y2_im.val[0]);
1075// INV_TOPJm(t_re.val[3], t_im.val[3], x2_re.val[1], x2_im.val[1], y2_re.val[1], y2_im.val[1]);
1076
1077// INV_BOTJ_LANE (y1_re.val[0], y1_im.val[0], t_re.val[0], t_im.val[0], s1_re_im);
1078// INV_BOTJ_LANE (y1_re.val[1], y1_im.val[1], t_re.val[1], t_im.val[1], s1_re_im);
1079
1080// INV_BOTJm_LANE(y2_re.val[0], y2_im.val[0], t_re.val[2], t_im.val[2], s1_re_im);
1081// INV_BOTJm_LANE(y2_re.val[1], y2_im.val[1], t_re.val[3], t_im.val[3], s1_re_im);
1082// /*
1083// * Level 7
1084// * ( 0, 64) - ( 32, 96)
1085// * ( 1, 65) - ( 33, 97)
1086// * ( 0, 64) + ( 32, 96)
1087// * ( 1, 65) + ( 33, 97)
1088// * ( 32, 96) = @ * ( 0, 1)
1089// * ( 33, 97) = @ * ( 0, 1)
1090// *
1091// * ( 2, 66) - ( 34, 98)
1092// * ( 3, 67) - ( 35, 99)
1093// * ( 2, 66) + ( 34, 98)
1094// * ( 3, 67) + ( 35, 99)
1095// * ( 34, 98) = @ * ( 0, 1)
1096// * ( 35, 99) = @ * ( 0, 1)
1097// * ----
1098// * ( 16, 80) - ( 48, 112)
1099// * ( 17, 81) - ( 49, 113)
1100// * ( 16, 80) + ( 48, 112)
1101// * ( 17, 81) + ( 49, 113)
1102// * ( 48, 112) = @ * ( 0, 1)
1103// * ( 49, 113) = @ * ( 0, 1)
1104// *
1105// * ( 18, 82) - ( 50, 114)
1106// * ( 19, 83) - ( 51, 115)
1107// * ( 18, 82) + ( 50, 114)
1108// * ( 19, 83) + ( 51, 115)
1109// * ( 50, 114) = @ * ( 0, 1)
1110// * ( 51, 115) = @ * ( 0, 1)
1111// */
1112
1113// if (Jm)
1114// {
1115// INV_TOPJm(t_re.val[0], t_im.val[0], x1_re.val[0], x1_im.val[0], x2_re.val[0], x2_im.val[0]);
1116// INV_TOPJm(t_re.val[1], t_im.val[1], x1_re.val[1], x1_im.val[1], x2_re.val[1], x2_im.val[1]);
1117
1118// INV_TOPJm(t_re.val[2], t_im.val[2], y1_re.val[0], y1_im.val[0], y2_re.val[0], y2_im.val[0]);
1119// INV_TOPJm(t_re.val[3], t_im.val[3], y1_re.val[1], y1_im.val[1], y2_re.val[1], y2_im.val[1]);
1120
1121// INV_BOTJm_LANE(x2_re.val[0], x2_im.val[0], t_re.val[0], t_im.val[0], s2_re_im);
1122// INV_BOTJm_LANE(x2_re.val[1], x2_im.val[1], t_re.val[1], t_im.val[1], s2_re_im);
1123// INV_BOTJm_LANE(y2_re.val[0], y2_im.val[0], t_re.val[2], t_im.val[2], s2_re_im);
1124// INV_BOTJm_LANE(y2_re.val[1], y2_im.val[1], t_re.val[3], t_im.val[3], s2_re_im);
1125// }
1126// else
1127// {
1128// INV_TOPJ(t_re.val[0], t_im.val[0], x1_re.val[0], x1_im.val[0], x2_re.val[0], x2_im.val[0]);
1129// INV_TOPJ(t_re.val[1], t_im.val[1], x1_re.val[1], x1_im.val[1], x2_re.val[1], x2_im.val[1]);
1130
1131// INV_TOPJ(t_re.val[2], t_im.val[2], y1_re.val[0], y1_im.val[0], y2_re.val[0], y2_im.val[0]);
1132// INV_TOPJ(t_re.val[3], t_im.val[3], y1_re.val[1], y1_im.val[1], y2_re.val[1], y2_im.val[1]);
1133
1134// INV_BOTJ_LANE(x2_re.val[0], x2_im.val[0], t_re.val[0], t_im.val[0], s2_re_im);
1135// INV_BOTJ_LANE(x2_re.val[1], x2_im.val[1], t_re.val[1], t_im.val[1], s2_re_im);
1136// INV_BOTJ_LANE(y2_re.val[0], y2_im.val[0], t_re.val[2], t_im.val[2], s2_re_im);
1137// INV_BOTJ_LANE(y2_re.val[1], y2_im.val[1], t_re.val[3], t_im.val[3], s2_re_im);
1138// }
1139
1140// vstorex2(&f[j + 2*len], x2_re);
1141// vstorex2(&f[j + 2*len + hn], x2_im);
1142
1143// vstorex2(&f[j + 3*len], y2_re);
1144// vstorex2(&f[j + 3*len + hn], y2_im);
1145
1146// if (!last)
1147// {
1148// vfmuln(x1_re.val[0], x1_re.val[0], fpr_p2_tab[logn]);
1149// vfmuln(x1_re.val[1], x1_re.val[1], fpr_p2_tab[logn]);
1150// vfmuln(x1_im.val[0], x1_im.val[0], fpr_p2_tab[logn]);
1151// vfmuln(x1_im.val[1], x1_im.val[1], fpr_p2_tab[logn]);
1152
1153// vfmuln(y1_re.val[0], y1_re.val[0], fpr_p2_tab[logn]);
1154// vfmuln(y1_re.val[1], y1_re.val[1], fpr_p2_tab[logn]);
1155// vfmuln(y1_im.val[0], y1_im.val[0], fpr_p2_tab[logn]);
1156// vfmuln(y1_im.val[1], y1_im.val[1], fpr_p2_tab[logn]);
1157// }
1158
1159// vstorex2(&f[j], x1_re);
1160// vstorex2(&f[j + hn], x1_im);
1161
1162// vstorex2(&f[j + len], y1_re);
1163// vstorex2(&f[j + len + hn], y1_im);
1164
1165// }
1166// }
1167// }
1168// }
1169
1170
1171
1172static
1173 void ZfN(iFFT_logn2)(fpr *f, const unsigned logn, const unsigned level, unsigned last)
1174{
1175 const unsigned int falcon_n = 1 << logn;
1176 const unsigned int hn = falcon_n >> 1;
1177
1178 // Total SIMD register: 26 = 16 + 8 + 2
1179 float64x2x4_t t_re, t_im; // 8
1180 float64x2x2_t x1_re, x2_re, x1_im, x2_im,
1181 y1_re, y2_re, y1_im, y2_im; // 16
1182 float64x2_t s1_re_im, s2_re_im; // 2
1183
1184 const fpr *fpr_inv_tab1 = NULL, *fpr_inv_tab2 = NULL;
1185 unsigned l, len, start, j, k1, k2;
1186 unsigned bar = logn - 4;
1187
1188 for (l = 4; l < logn - level - 1; l += 2)
1189 {
1190 len = 1 << l;
1191 last -= 1;
1192 fpr_inv_tab1 = fpr_table[bar--];
1193 fpr_inv_tab2 = fpr_table[bar--];
1194 k1 = 0; k2 = 0;
1195
1196 for (start = 0; start < hn; start += 1 << (l + 2))
1197 {
1198 vload(s1_re_im, &fpr_inv_tab1[k1]);
1199 vload(s2_re_im, &fpr_inv_tab2[k2]);
1200 k1 += 2;
1201 k2 += 2 * ((start & 127) == 64);
1202 if (!last)
1203 {
1204 vfmuln(s2_re_im, s2_re_im, fpr_p2_tab[logn]);
1205 }
1206 for (j = start; j < start + len; j += 4)
1207 {
1208 /*
1209 Level 6
1210 * ( 0, 64) - ( 16, 80)
1211 * ( 1, 65) - ( 17, 81)
1212 * ( 0, 64) + ( 16, 80)
1213 * ( 1, 65) + ( 17, 81)
1214 * ( 16, 80) = @ * ( 0, 1)
1215 * ( 17, 81) = @ * ( 0, 1)
1216 *
1217 * ( 2, 66) - ( 18, 82)
1218 * ( 3, 67) - ( 19, 83)
1219 * ( 2, 66) + ( 18, 82)
1220 * ( 3, 67) + ( 19, 83)
1221 * ( 18, 82) = @ * ( 0, 1)
1222 * ( 19, 83) = @ * ( 0, 1)
1223 *
1224 * ( 32, 96) - ( 48, 112)
1225 * ( 33, 97) - ( 49, 113)
1226 * ( 32, 96) + ( 48, 112)
1227 * ( 33, 97) + ( 49, 113)
1228 * ( 48, 112) = j@ * ( 0, 1)
1229 * ( 49, 113) = j@ * ( 0, 1)
1230 *
1231 * ( 34, 98) - ( 50, 114)
1232 * ( 35, 99) - ( 51, 115)
1233 * ( 34, 98) + ( 50, 114)
1234 * ( 35, 99) + ( 51, 115)
1235 * ( 50, 114) = j@ * ( 0, 1)
1236 * ( 51, 115) = j@ * ( 0, 1)
1237 */
1238 // x1: 0 -> 4 | 64 -> 67
1239 // y1: 16 -> 19 | 80 -> 81
1240 // x2: 32 -> 35 | 96 -> 99
1241 // y2: 48 -> 51 | 112 -> 115
1242 vloadx2(x1_re, &f[j]);
1243 vloadx2(x1_im, &f[j + hn]);
1244 vloadx2(y1_re, &f[j + len]);
1245 vloadx2(y1_im, &f[j + len + hn]);
1246
1247 INV_TOPJ (t_re.val[0], t_im.val[0], x1_re.val[0], x1_im.val[0], y1_re.val[0], y1_im.val[0]);
1248 INV_TOPJ (t_re.val[1], t_im.val[1], x1_re.val[1], x1_im.val[1], y1_re.val[1], y1_im.val[1]);
1249
1250 vloadx2(x2_re, &f[j + 2*len]);
1251 vloadx2(x2_im, &f[j + 2*len + hn]);
1252 vloadx2(y2_re, &f[j + 3*len]);
1253 vloadx2(y2_im, &f[j + 3*len + hn]);
1254
1255 INV_TOPJm(t_re.val[2], t_im.val[2], x2_re.val[0], x2_im.val[0], y2_re.val[0], y2_im.val[0]);
1256 INV_TOPJm(t_re.val[3], t_im.val[3], x2_re.val[1], x2_im.val[1], y2_re.val[1], y2_im.val[1]);
1257
1258 INV_BOTJ_LANE (y1_re.val[0], y1_im.val[0], t_re.val[0], t_im.val[0], s1_re_im);
1259 INV_BOTJ_LANE (y1_re.val[1], y1_im.val[1], t_re.val[1], t_im.val[1], s1_re_im);
1260
1261 INV_BOTJm_LANE(y2_re.val[0], y2_im.val[0], t_re.val[2], t_im.val[2], s1_re_im);
1262 INV_BOTJm_LANE(y2_re.val[1], y2_im.val[1], t_re.val[3], t_im.val[3], s1_re_im);
1263 /*
1264 * Level 7
1265 * ( 0, 64) - ( 32, 96)
1266 * ( 1, 65) - ( 33, 97)
1267 * ( 0, 64) + ( 32, 96)
1268 * ( 1, 65) + ( 33, 97)
1269 * ( 32, 96) = @ * ( 0, 1)
1270 * ( 33, 97) = @ * ( 0, 1)
1271 *
1272 * ( 2, 66) - ( 34, 98)
1273 * ( 3, 67) - ( 35, 99)
1274 * ( 2, 66) + ( 34, 98)
1275 * ( 3, 67) + ( 35, 99)
1276 * ( 34, 98) = @ * ( 0, 1)
1277 * ( 35, 99) = @ * ( 0, 1)
1278 * ----
1279 * ( 16, 80) - ( 48, 112)
1280 * ( 17, 81) - ( 49, 113)
1281 * ( 16, 80) + ( 48, 112)
1282 * ( 17, 81) + ( 49, 113)
1283 * ( 48, 112) = @ * ( 0, 1)
1284 * ( 49, 113) = @ * ( 0, 1)
1285 *
1286 * ( 18, 82) - ( 50, 114)
1287 * ( 19, 83) - ( 51, 115)
1288 * ( 18, 82) + ( 50, 114)
1289 * ( 19, 83) + ( 51, 115)
1290 * ( 50, 114) = @ * ( 0, 1)
1291 * ( 51, 115) = @ * ( 0, 1)
1292 */
1293
1294
1295 INV_TOPJ(t_re.val[0], t_im.val[0], x1_re.val[0], x1_im.val[0], x2_re.val[0], x2_im.val[0]);
1296 INV_TOPJ(t_re.val[1], t_im.val[1], x1_re.val[1], x1_im.val[1], x2_re.val[1], x2_im.val[1]);
1297
1298 INV_TOPJ(t_re.val[2], t_im.val[2], y1_re.val[0], y1_im.val[0], y2_re.val[0], y2_im.val[0]);
1299 INV_TOPJ(t_re.val[3], t_im.val[3], y1_re.val[1], y1_im.val[1], y2_re.val[1], y2_im.val[1]);
1300
1301 INV_BOTJ_LANE(x2_re.val[0], x2_im.val[0], t_re.val[0], t_im.val[0], s2_re_im);
1302 INV_BOTJ_LANE(x2_re.val[1], x2_im.val[1], t_re.val[1], t_im.val[1], s2_re_im);
1303 INV_BOTJ_LANE(y2_re.val[0], y2_im.val[0], t_re.val[2], t_im.val[2], s2_re_im);
1304 INV_BOTJ_LANE(y2_re.val[1], y2_im.val[1], t_re.val[3], t_im.val[3], s2_re_im);
1305
1306 vstorex2(&f[j + 2*len], x2_re);
1307 vstorex2(&f[j + 2*len + hn], x2_im);
1308
1309 vstorex2(&f[j + 3*len], y2_re);
1310 vstorex2(&f[j + 3*len + hn], y2_im);
1311
1312 if (!last)
1313 {
1314 vfmuln(x1_re.val[0], x1_re.val[0], fpr_p2_tab[logn]);
1315 vfmuln(x1_re.val[1], x1_re.val[1], fpr_p2_tab[logn]);
1316 vfmuln(x1_im.val[0], x1_im.val[0], fpr_p2_tab[logn]);
1317 vfmuln(x1_im.val[1], x1_im.val[1], fpr_p2_tab[logn]);
1318
1319 vfmuln(y1_re.val[0], y1_re.val[0], fpr_p2_tab[logn]);
1320 vfmuln(y1_re.val[1], y1_re.val[1], fpr_p2_tab[logn]);
1321 vfmuln(y1_im.val[0], y1_im.val[0], fpr_p2_tab[logn]);
1322 vfmuln(y1_im.val[1], y1_im.val[1], fpr_p2_tab[logn]);
1323 }
1324
1325 vstorex2(&f[j], x1_re);
1326 vstorex2(&f[j + hn], x1_im);
1327
1328 vstorex2(&f[j + len], y1_re);
1329 vstorex2(&f[j + len + hn], y1_im);
1330
1331 }
1332 //
1333 start += 1 << (l + 2);
1334 if (start >= hn) break;
1335
1336 vload(s1_re_im, &fpr_inv_tab1[k1]);
1337 vload(s2_re_im, &fpr_inv_tab2[k2]);
1338 k1 += 2;
1339 k2 += 2 * ((start & 127) == 64);
1340 if (!last)
1341 {
1342 vfmuln(s2_re_im, s2_re_im, fpr_p2_tab[logn]);
1343 }
1344
1345 for (j = start; j < start + len; j += 4)
1346 {
1347 /*
1348 Level 6
1349 * ( 0, 64) - ( 16, 80)
1350 * ( 1, 65) - ( 17, 81)
1351 * ( 0, 64) + ( 16, 80)
1352 * ( 1, 65) + ( 17, 81)
1353 * ( 16, 80) = @ * ( 0, 1)
1354 * ( 17, 81) = @ * ( 0, 1)
1355 *
1356 * ( 2, 66) - ( 18, 82)
1357 * ( 3, 67) - ( 19, 83)
1358 * ( 2, 66) + ( 18, 82)
1359 * ( 3, 67) + ( 19, 83)
1360 * ( 18, 82) = @ * ( 0, 1)
1361 * ( 19, 83) = @ * ( 0, 1)
1362 *
1363 * ( 32, 96) - ( 48, 112)
1364 * ( 33, 97) - ( 49, 113)
1365 * ( 32, 96) + ( 48, 112)
1366 * ( 33, 97) + ( 49, 113)
1367 * ( 48, 112) = j@ * ( 0, 1)
1368 * ( 49, 113) = j@ * ( 0, 1)
1369 *
1370 * ( 34, 98) - ( 50, 114)
1371 * ( 35, 99) - ( 51, 115)
1372 * ( 34, 98) + ( 50, 114)
1373 * ( 35, 99) + ( 51, 115)
1374 * ( 50, 114) = j@ * ( 0, 1)
1375 * ( 51, 115) = j@ * ( 0, 1)
1376 */
1377 // x1: 0 -> 4 | 64 -> 67
1378 // y1: 16 -> 19 | 80 -> 81
1379 // x2: 32 -> 35 | 96 -> 99
1380 // y2: 48 -> 51 | 112 -> 115
1381 vloadx2(x1_re, &f[j]);
1382 vloadx2(x1_im, &f[j + hn]);
1383 vloadx2(y1_re, &f[j + len]);
1384 vloadx2(y1_im, &f[j + len + hn]);
1385
1386 INV_TOPJ (t_re.val[0], t_im.val[0], x1_re.val[0], x1_im.val[0], y1_re.val[0], y1_im.val[0]);
1387 INV_TOPJ (t_re.val[1], t_im.val[1], x1_re.val[1], x1_im.val[1], y1_re.val[1], y1_im.val[1]);
1388
1389 vloadx2(x2_re, &f[j + 2*len]);
1390 vloadx2(x2_im, &f[j + 2*len + hn]);
1391 vloadx2(y2_re, &f[j + 3*len]);
1392 vloadx2(y2_im, &f[j + 3*len + hn]);
1393
1394 INV_TOPJm(t_re.val[2], t_im.val[2], x2_re.val[0], x2_im.val[0], y2_re.val[0], y2_im.val[0]);
1395 INV_TOPJm(t_re.val[3], t_im.val[3], x2_re.val[1], x2_im.val[1], y2_re.val[1], y2_im.val[1]);
1396
1397 INV_BOTJ_LANE (y1_re.val[0], y1_im.val[0], t_re.val[0], t_im.val[0], s1_re_im);
1398 INV_BOTJ_LANE (y1_re.val[1], y1_im.val[1], t_re.val[1], t_im.val[1], s1_re_im);
1399
1400 INV_BOTJm_LANE(y2_re.val[0], y2_im.val[0], t_re.val[2], t_im.val[2], s1_re_im);
1401 INV_BOTJm_LANE(y2_re.val[1], y2_im.val[1], t_re.val[3], t_im.val[3], s1_re_im);
1402 /*
1403 * Level 7
1404 * ( 0, 64) - ( 32, 96)
1405 * ( 1, 65) - ( 33, 97)
1406 * ( 0, 64) + ( 32, 96)
1407 * ( 1, 65) + ( 33, 97)
1408 * ( 32, 96) = @ * ( 0, 1)
1409 * ( 33, 97) = @ * ( 0, 1)
1410 *
1411 * ( 2, 66) - ( 34, 98)
1412 * ( 3, 67) - ( 35, 99)
1413 * ( 2, 66) + ( 34, 98)
1414 * ( 3, 67) + ( 35, 99)
1415 * ( 34, 98) = @ * ( 0, 1)
1416 * ( 35, 99) = @ * ( 0, 1)
1417 * ----
1418 * ( 16, 80) - ( 48, 112)
1419 * ( 17, 81) - ( 49, 113)
1420 * ( 16, 80) + ( 48, 112)
1421 * ( 17, 81) + ( 49, 113)
1422 * ( 48, 112) = @ * ( 0, 1)
1423 * ( 49, 113) = @ * ( 0, 1)
1424 *
1425 * ( 18, 82) - ( 50, 114)
1426 * ( 19, 83) - ( 51, 115)
1427 * ( 18, 82) + ( 50, 114)
1428 * ( 19, 83) + ( 51, 115)
1429 * ( 50, 114) = @ * ( 0, 1)
1430 * ( 51, 115) = @ * ( 0, 1)
1431 */
1432
1433 INV_TOPJm(t_re.val[0], t_im.val[0], x1_re.val[0], x1_im.val[0], x2_re.val[0], x2_im.val[0]);
1434 INV_TOPJm(t_re.val[1], t_im.val[1], x1_re.val[1], x1_im.val[1], x2_re.val[1], x2_im.val[1]);
1435
1436 INV_TOPJm(t_re.val[2], t_im.val[2], y1_re.val[0], y1_im.val[0], y2_re.val[0], y2_im.val[0]);
1437 INV_TOPJm(t_re.val[3], t_im.val[3], y1_re.val[1], y1_im.val[1], y2_re.val[1], y2_im.val[1]);
1438
1439 INV_BOTJm_LANE(x2_re.val[0], x2_im.val[0], t_re.val[0], t_im.val[0], s2_re_im);
1440 INV_BOTJm_LANE(x2_re.val[1], x2_im.val[1], t_re.val[1], t_im.val[1], s2_re_im);
1441 INV_BOTJm_LANE(y2_re.val[0], y2_im.val[0], t_re.val[2], t_im.val[2], s2_re_im);
1442 INV_BOTJm_LANE(y2_re.val[1], y2_im.val[1], t_re.val[3], t_im.val[3], s2_re_im);
1443
1444 vstorex2(&f[j + 2*len], x2_re);
1445 vstorex2(&f[j + 2*len + hn], x2_im);
1446
1447 vstorex2(&f[j + 3*len], y2_re);
1448 vstorex2(&f[j + 3*len + hn], y2_im);
1449
1450 if (!last)
1451 {
1452 vfmuln(x1_re.val[0], x1_re.val[0], fpr_p2_tab[logn]);
1453 vfmuln(x1_re.val[1], x1_re.val[1], fpr_p2_tab[logn]);
1454 vfmuln(x1_im.val[0], x1_im.val[0], fpr_p2_tab[logn]);
1455 vfmuln(x1_im.val[1], x1_im.val[1], fpr_p2_tab[logn]);
1456
1457 vfmuln(y1_re.val[0], y1_re.val[0], fpr_p2_tab[logn]);
1458 vfmuln(y1_re.val[1], y1_re.val[1], fpr_p2_tab[logn]);
1459 vfmuln(y1_im.val[0], y1_im.val[0], fpr_p2_tab[logn]);
1460 vfmuln(y1_im.val[1], y1_im.val[1], fpr_p2_tab[logn]);
1461 }
1462
1463 vstorex2(&f[j], x1_re);
1464 vstorex2(&f[j + hn], x1_im);
1465
1466 vstorex2(&f[j + len], y1_re);
1467 vstorex2(&f[j + len + hn], y1_im);
1468
1469 }
1470 //
1471 }
1472 }
1473}
1474
1475
1476/*
1477* Support logn from [1, 10]
1478* Can be easily extended to logn > 10
1479*/
1480void ZfN(iFFT)(fpr *f, const unsigned logn)
1481{
1482 const unsigned level = (logn - 5) & 1;
1483
1484 switch (logn)
1485 {
1486 case 2:
1487 ZfN(iFFT_log2)(f);
1488 break;
1489
1490 case 3:
1491 ZfN(iFFT_log3)(f);
1492 break;
1493
1494 case 4:
1495 ZfN(iFFT_log4)(f);
1496 break;
1497
1498 case 5:
1499 ZfN(iFFT_log5)(f, 5, 1);
1500 break;
1501
1502 case 6:
1503 ZfN(iFFT_log5)(f, logn, 0);
1504 ZfN(iFFT_logn1)(f, logn, 1);
1505 break;
1506
1507 case 7:
1508 case 9:
1509 ZfN(iFFT_log5)(f, logn, 0);
1510 ZfN(iFFT_logn2)(f, logn, level, 1);
1511 break;
1512
1513 case 8:
1514 case 10:
1515 ZfN(iFFT_log5)(f, logn, 0);
1516 ZfN(iFFT_logn2)(f, logn, level, 0);
1517 ZfN(iFFT_logn1)(f, logn, 1);
1518 break;
1519
1520 default:
1521 break;
1522 }
1523}
1524#endif
1525
1526// generic fft stuff
1527
1528#include "../commons_private.h"
1529#include "reim_fft_internal.h"
1530#include "reim_fft_private.h"
1531
1532void reim_ctwiddle(double* ra, double* ia, double* rb, double* ib, double omre, double omim);
1533void reim_citwiddle(double* ra, double* ia, double* rb, double* ib, double omre, double omim);
1534
1535void reim_fft16_ref(double* dre, double* dim, const void* pom);
1536void fill_reim_fft16_omegas(const double entry_pwr, double** omg);
1537void reim_fft8_ref(double* dre, double* dim, const void* pom);
1538void reim_twiddle_fft_ref(uint64_t h, double* re, double* im, double om[2]);
1539void fill_reim_twiddle_fft_ref(const double s, double** omg);
1540void reim_bitwiddle_fft_ref(uint64_t h, double* re, double* im, double om[4]);
1541void fill_reim_bitwiddle_fft_ref(const double s, double** omg);
1542void reim_fft_rec_16_ref(uint64_t m, double* re, double* im, double** omg);
1543void fill_reim_fft_rec_16_omegas(uint64_t m, double entry_pwr, double** omg);
1544
1545void fill_reim_twiddle_fft_omegas_ref(const double rs0, double** omg) {
1546 (*omg)[0] = cos(2 * M_PI * rs0);
1547 (*omg)[1] = sin(2 * M_PI * rs0);
1548 *omg += 2;
1549}
1550
1551void fill_reim_bitwiddle_fft_omegas_ref(const double rs0, double** omg) {
1552 double rs1 = 2. * rs0;
1553 (*omg)[0] = cos(2 * M_PI * rs1);
1554 (*omg)[1] = sin(2 * M_PI * rs1);
1555 (*omg)[2] = cos(2 * M_PI * rs0);
1556 (*omg)[3] = sin(2 * M_PI * rs0);
1557 *omg += 4;
1558}
1559
1560#define reim_fft16_f reim_fft16_neon
1561#define reim_fft16_pom_offset 16
1562#define fill_reim_fft16_omegas_f fill_reim_fft16_omegas_neon
1563// currently, m=4 uses the ref implem (corner case with low impact) TODO!!
1564#define reim_fft8_f reim_fft8_ref
1565#define reim_fft8_pom_offset 8
1566#define fill_reim_fft8_omegas_f fill_reim_fft8_omegas
1567// currently, m=4 uses the ref implem (corner case with low impact) TODO!!
1568#define reim_fft4_f reim_fft4_ref
1569#define reim_fft4_pom_offset 4
1570#define fill_reim_fft4_omegas_f fill_reim_fft4_omegas
1571// m = 2 will use the ref implem, since intrinsics don't provide any speed-up
1572#define reim_fft2_f reim_fft2_ref
1573#define reim_fft2_pom_offset 2
1574#define fill_reim_fft2_omegas_f fill_reim_fft2_omegas
1575
1576// neon twiddle use the same omegas layout as the ref implem
1577#define reim_twiddle_fft_f reim_twiddle_fft_neon
1578#define reim_twiddle_fft_pom_offset 2
1579#define fill_reim_twiddle_fft_omegas_f fill_reim_twiddle_fft_omegas_ref
1580
1581// neon bi-twiddle use the same omegas layout as the ref implem
1582#define reim_bitwiddle_fft_f reim_bitwiddle_fft_neon
1583#define reim_bitwiddle_fft_pom_offset 4
1584#define fill_reim_bitwiddle_fft_omegas_f fill_reim_bitwiddle_fft_omegas_ref
1585
1586// template functions to produce
1587#define reim_fft_bfs_16_f reim_fft_bfs_16_neon
1588#define fill_reim_fft_bfs_16_omegas_f fill_reim_fft_bfs_16_omegas_neon
1589#define reim_fft_rec_16_f reim_fft_rec_16_neon
1590#define fill_reim_fft_rec_16_omegas_f fill_reim_fft_rec_16_omegas_neon
1591#define reim_fft_f reim_fft_neon
1592#define fill_reim_fft_omegas_f fill_reim_fft_omegas_neon
1593
1594#include "reim_fft_core_template.h"
1595
1596EXPORT REIM_FFT_PRECOMP* new_reim_fft_precomp_neon(uint32_t m, uint32_t num_buffers) {
1597 const uint64_t OMG_SPACE = ceilto64b(2 * m * sizeof(double));
1598 const uint64_t BUF_SIZE = ceilto64b(2 * m * sizeof(double));
1599 void* reps = malloc(sizeof(REIM_FFT_PRECOMP) // base
1600 + 63 // padding
1601 + OMG_SPACE // tables //TODO 16?
1602 + num_buffers * BUF_SIZE // buffers
1603 );
1604 uint64_t aligned_addr = ceilto64b((uint64_t)(reps) + sizeof(REIM_FFT_PRECOMP));
1605 REIM_FFT_PRECOMP* r = (REIM_FFT_PRECOMP*)reps;
1606 r->m = m;
1607 r->buf_size = BUF_SIZE;
1608 r->powomegas = (double*)aligned_addr;
1609 r->aligned_buffers = (void*)(aligned_addr + OMG_SPACE);
1610 // fill in powomegas
1611 double* omg = (double*)r->powomegas;
1612 fill_reim_fft_omegas_f(m, 0.25, &omg);
1613 if (((uint64_t)omg) - aligned_addr > OMG_SPACE) abort();
1614 // dispatch the right implementation
1615 //{
1616 // if (CPU_SUPPORTS("fma")) {
1617 // r->function = reim_fft_avx2_fma;
1618 // } else {
1619 // r->function = reim_fft_ref;
1620 // }
1621 //}
1622 r->function = reim_fft_f;
1623 return reps;
1624}