1 /*****************************************************************************
2 The functions in this file are adaptations of csfft.c that should be
3 thread-safe for use with OpenMP. As such, this file is intended to be
4 #include-d into a source file that uses OpenMP directives. In such a case,
5 the following inclusions are not really necessary, but included here for
6 reference.
7
8 When using these functions, csfft_cox_OMP_SETUP() must be called once
9 first, in order to setup some thread specific static storage. These
10 variable names all start with 'CFO_'.
11 -- RWCox -- Dec 2015
12
13 ***** I'm not actually sure these functions work -- BEWARE !! :( *****
14
15 ***** Instead ... try using the contents of file fftn_OMP.c *****
16 ******************************************************************************/
17
18 #include "mrilib.h" /* AFNI package library header */
19 #ifdef USE_OMP
20 # include <omp.h>
21 #endif
22
23 /***=====================================================================**
24 *** Prototypes for externally callable functions: **
25 *** Complex-to-complex FFT in place: **
26 *** mode = -1 or +1 (for inverse scaling, csfft_scale_inverse_OMP **
27 *** idim = dimension (power of 2, maybe with factors of 3 and 5) **
28 *** xc = input/output array **
29 *** Re-initializes itself only when idim changes from previous call. **
30 ***=====================================================================**/
31
32 void csfft_cox_OMP( int mode , int idim , complex *xc ) ;
33 void csfft_scale_inverse_OMP( int scl ) ; /* scl=1 ==> force 1/N for mode=+1 **/
34 /* scl=0 ==> no 1/N scaling **/
35
36 /*** Aug 1999: **
37 *** idim can now contain factors of 3 and/or 5, up to and including **
38 *** 3^3 * 5^3. Some examples: **
39 *** 135 144 150 160 180 192 200 216 225 240 250 **
40 *** The routine csfft_nextup(n) returns the smallest size FFT **
41 *** >= n that csfft_cox_OMP() knows how to do. **
42 *** Note that efficiency of the different lengths can be quite **
43 *** different. In general, powers of 2 are fastest, but a **
44 *** single factor of 3 and/or 5 doesn't cause too much slowdown. **
45 *** The routine csfft_nextup_one35(n) returns the smallest size FFT **
46 *** >= n that contains at most one power of 3, at most one power **
47 *** of 5, and least one power of 2. In trials, these are the most **
48 *** time efficient. **/
49
50 /*-- Aug 1999: internal routines to do FFTs by decimation by 3, 4, 5 --*/
51
52 static void fft_3dec_OMP( int , int , complex * ) ;
53 static void fft_4dec_OMP( int , int , complex * ) ;
54 static void fft_5dec_OMP( int , int , complex * ) ;
55
56 #undef RMAX
57 #define RMAX 7 /* maximum levels of recursion for decimation by 3, 4, 5 */
58
59 #undef PI
60 #define PI (3.141592653589793238462643) /* that should be good enough */
61
62 /*---------------------------------------------------------------------------*/
63 /* When using these OpenMP-ized functions,
64 csfft_cox_OMP_SETUP() must be called once first! [Dec 2015]
65 *//*-------------------------------------------------------------------------*/
66
67 static int CFO_nthr = 1 ; /* max number of threads to allow for */
68 static complex **CFO_csplus =NULL ; /* per-thread arrays of trig consts */
69 static complex **CFO_csminus=NULL ;
70 static int *CFO_nold =NULL ; /* per-thread last dimension used */
71 static int *CFO_sclinv =NULL ; /* per-thread inverse scaling flag */
72 static int *CFO_recC =NULL ; /* per-thread recursion level for csfft_cox_OMP */
73 static int *CFO_recX =NULL ; /* per-thread recursion level for decimaion funcs */
74
75 static int *CFO_rmold3[RMAX] ; /* per-thread, per-recursion last M value used */
76 static int *CFO_rmold4[RMAX] ; /* for decimation by 3, 4, 5 functions */
77 static int *CFO_rmold5[RMAX] ;
78
79 static int *nCFO_aa[RMAX] ; /* per-thread, per-recursion temp arrays */
80 static int *nCFO_bb[RMAX] ; /* nCFO_aa[r][i] is the size of the 'aa' */
81 static int *nCFO_cc[RMAX] ; /* array for recursion level r, thread i */
82 static int *nCFO_dd[RMAX] ; /* CFO_aa[r][i] is the corresponding array */
83 static int *nCFO_ee[RMAX] ;
84 static int *nCFO_cs[RMAX] ;
85 static complex **CFO_aa[RMAX] ;
86 static complex **CFO_bb[RMAX] ;
87 static complex **CFO_cc[RMAX] ;
88 static complex **CFO_dd[RMAX] ;
89 static complex **CFO_ee[RMAX] ;
90 static complex **CFO_cs[RMAX] ;
91
92 /* macro to declare and set the ithr = thread number (local variable) */
93
94 #undef DECLARE_ithr
95 #ifdef USE_OMP
96 # define DECLARE_ithr int ithr=omp_get_thread_num()
97 #else
98 # define DECLARE_ithr int ithr=0
99 #endif
100
101 /*----------------------------------------------------------------------------*/
102 /* Allocate the CFO arrays, just one time. */
103
csfft_cox_OMP_SETUP(void)104 void csfft_cox_OMP_SETUP(void)
105 {
106 int rr ;
107
108 if( CFO_nold != NULL ) return ; /* was called before! */
109
110 #ifdef USE_OMP
111 CFO_nthr = omp_get_max_threads() ;
112 #endif
113
114 /* arrays that don't depend on recursion level, just thread number */
115
116 CFO_nold = (int *) calloc(sizeof(int) ,CFO_nthr) ;
117 CFO_sclinv = (int *) calloc(sizeof(int) ,CFO_nthr) ;
118 CFO_recC = (int *) calloc(sizeof(int) ,CFO_nthr) ;
119 CFO_recX = (int *) calloc(sizeof(int) ,CFO_nthr) ;
120 CFO_csplus = (complex **)calloc(sizeof(complex *),CFO_nthr) ;
121 CFO_csminus = (complex **)calloc(sizeof(complex *),CFO_nthr) ;
122
123 /* arrays that depend on recursion level and on thread number;
124 note that we are just allocating pointer space now; the
125 actual arrays are allocated by the CFO_assign macro (infra) */
126
127 for( rr=0 ; rr < RMAX ; rr++ ){
128 nCFO_aa[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ;
129 nCFO_bb[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ;
130 nCFO_cc[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ;
131 nCFO_dd[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ;
132 nCFO_ee[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ;
133 nCFO_cs[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ;
134 CFO_aa[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ;
135 CFO_bb[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ;
136 CFO_cc[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ;
137 CFO_dd[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ;
138 CFO_ee[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ;
139 CFO_cs[rr] = (complex **)calloc(sizeof(complex *),CFO_nthr) ;
140
141 CFO_rmold3[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ;
142 CFO_rmold4[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ;
143 CFO_rmold5[rr] = (int *) calloc(sizeof(int) ,CFO_nthr) ;
144 }
145
146 return ;
147 }
148
149 /** this macro assigns a workspace of length ll (complex), at
150 recursion level rec (0..RMAX-1), to pointer named 'ab', from
151 static array CFO_'ab'[rr][ithr], in thread number ithr.
152 ithr and rec are local to the function whence this is called!
153 The name 'ab' can be one of 'aa' 'bb' 'cc' 'dd' 'ee' 'cs'.
154 The C preprocessor token catenation operator ## is used here. */
155
156 #undef CFO_assign
157 #define CFO_assign(ab,ll) \
158 do{ if( nCFO_##ab[rec][ithr] < (ll) ){ \
159 CFO_##ab[rec][ithr] = (complex *)realloc(CFO_##ab[rec][ithr], \
160 sizeof(complex)*(ll)); \
161 nCFO_##ab[rec][ithr] = (ll) ; \
162 } \
163 ab = CFO_##ab[rec][ithr] ; \
164 } while(0)
165
166 /*-------------- To order the 1/N scaling on inversion: Aug 1999 ------------*/
167
csfft_scale_inverse_OMP(int scl)168 void csfft_scale_inverse_OMP( int scl ){
169 DECLARE_ithr ; CFO_sclinv[ithr] = scl; return;
170 }
171
172 /*-------------- For the unrolled FFT routines: November 1998 --------------*/
173
174 static void fft8 ( int mode , complex *xc ) ; /* completely */
175 static void fft16( int mode , complex *xc ) ; /* unrolled */
176 static void fft32( int mode , complex *xc ) ; /* routines */
177 static void fft64( int mode , complex *xc ) ;
178
179 /* convenience macros for decimation by 4 */
180
181 #undef DONT_USE_4DEC
182 #ifndef DONT_USE_4DEC
183 #define fft128(m,x) fft_4dec_OMP(m, 128,x) /* via fft32 (non-recursive) */
184 #define fft256(m,x) fft_4dec_OMP(m, 256,x) /* via fft64 (non-recursive) */
185 #define fft512(m,x) fft_4dec_OMP(m, 512,x) /* via fft128 (recursive X 1) */
186 #define fft1024(m,x) fft_4dec_OMP(m, 1024,x) /* via fft256 (recursive X 1) */
187 #if 1
188 #define fft2048(m,x) fft_4dec_OMP(m, 2048,x) /* via fft512 (recursive X 2) */
189 #define fft4096(m,x) fft_4dec_OMP(m, 4096,x) /* via fft1024 (recursive X 2) */
190 #define fft8192(m,x) fft_4dec_OMP(m, 8192,x) /* via fft2048 (recursive X 3) */
191 #define fft16384(m,x) fft_4dec_OMP(m,16384,x) /* via fft4096 (recursive X 3) */
192 #define fft32768(m,x) fft_4dec_OMP(m,32768,x) /* via fft8192 (recursive X 4) */
193 #define fft65536(m,x) fft_4dec_OMP(m,65536,x) /* via fft16384 (recursive X 4) */
194 #endif
195 #else
196 #define fft128(m,x) csfft_cox_OMP(m, 128,x)
197 #define fft256(m,x) csfft_cox_OMP(m, 256,x)
198 #define fft512(m,x) csfft_cox_OMP(m, 512,x)
199 #define fft1024(m,x) csfft_cox_OMP(m, 1024,x)
200 #if 1
201 #define fft2048(m,x) csfft_cox_OMP(m, 2048,x)
202 #define fft4096(m,x) csfft_cox_OMP(m, 4096,x)
203 #define fft8192(m,x) csfft_cox_OMP(m, 8192,x)
204 #define fft16384(m,x) csfft_cox_OMP(m,16384,x)
205 #define fft32768(m,x) csfft_cox_OMP(m,32768,x)
206 #define fft65536(m,x) csfft_cox_OMP(m,65536,x)
207 #endif
208 #endif
209
210 /*-- do FFTs of size 2 and 4 with macros, not functions --*/
211
212 #define fft2(m,x) do{ float a,b,c,d ; \
213 a = x[0].r + x[1].r ; b = x[0].i + x[1].i ; \
214 c = x[0].r - x[1].r ; d = x[0].i - x[1].i ; \
215 x[0].r = a ; x[0].i = b ; \
216 x[1].r = c ; x[1].i = d ; } while(0)
217
218 #define fft4(m,x) do{ float acpr,acmr,bdpr,bdmr, acpi,acmi,bdpi,bdmi; \
219 acpr = x[0].r + x[2].r; acmr = x[0].r - x[2].r; \
220 bdpr = x[1].r + x[3].r; bdmr = x[1].r - x[3].r; \
221 acpi = x[0].i + x[2].i; acmi = x[0].i - x[2].i; \
222 bdpi = x[1].i + x[3].i; bdmi = x[1].i - x[3].i; \
223 x[0].r = acpr+bdpr ; x[0].i = acpi+bdpi ; \
224 x[2].r = acpr-bdpr ; x[2].i = acpi-bdpi ; \
225 if(m > 0){ \
226 x[1].r = acmr-bdmi ; x[1].i = acmi+bdmr ; \
227 x[3].r = acmr+bdmi ; x[3].i = acmi-bdmr ; \
228 } else { \
229 x[1].r = acmr+bdmi ; x[1].i = acmi-bdmr ; \
230 x[3].r = acmr-bdmi ; x[3].i = acmi+bdmr ; \
231 } } while(0)
232
233 /*--------------------------------------------------------------------
234 Initialize csfft trig constants table. Adapted from AJ's code.
235 Does both mode=1 and mode=-1 tables. (mode=+1 == inversion)
236 ----------------------------------------------------------------------*/
237
csfft_trigconsts(int idim)238 static void csfft_trigconsts( int idim ) /* internal function */
239 {
240 DECLARE_ithr ;
241 register unsigned int m, n, i0, i1, i2, i3, k;
242 register complex *r0, *csp;
243 register float co, si, f0, f1, f2, f3, f4;
244 double al;
245 complex *my_csp , *my_csm ;
246
247 if( idim == CFO_nold[ithr] ) return ; /* previously done */
248
249 { if( idim > CFO_nold[ithr] ){
250 CFO_csplus[ithr] = (complex *)realloc( CFO_csplus[ithr] , sizeof(complex)*idim ) ;
251 CFO_csminus[ithr] = (complex *)realloc( CFO_csminus[ithr], sizeof(complex)*idim ) ;
252 if( CFO_csplus[ithr] == NULL || CFO_csminus[ithr] == NULL ){
253 fprintf(stderr,"\n*** csfft_cox_OMP cannot realloc space idim=%d! ***\n",idim); EXIT(1) ;
254 }
255 }}
256
257 CFO_nold[ithr] = n = idim ;
258
259 my_csp = CFO_csplus[ithr] ; my_csm = CFO_csminus[ithr] ;
260
261 f1 = 1.0 ; /* csplus init */
262 m = 1; k = 0;
263 while (n > m) {
264 i3 = m << 1; f2 = m; al = f1*PI/f2;
265 co = cos(al); si = sin(al);
266 (my_csp + k)->r = 1.; (my_csp + k)->i = 0.;
267 for (i0=0; i0 < m; i0++) {
268 k++;
269 csp = my_csp + k; r0 = csp - 1;
270 csp->r = r0->r * co - r0->i * si;
271 csp->i = r0->i * co + r0->r * si;
272 }
273 m = i3;
274 }
275
276 f1 = -1.0 ; /* csminus init */
277 m = 1; k = 0;
278 while (n > m) {
279 i3 = m << 1; f2 = m; al = f1*PI/f2;
280 co = cos(al); si = sin(al);
281 (my_csm + k)->r = 1.; (my_csm + k)->i = 0.;
282 for (i0=0; i0 < m; i0++) {
283 k++;
284 csp = my_csm + k; r0 = csp - 1;
285 csp->r = r0->r * co - r0->i * si;
286 csp->i = r0->i * co + r0->r * si;
287 }
288 m = i3;
289 }
290 return ;
291 }
292
293 /*--------------------------------------------------------------------
294 Complex-to-complex FFT in place:
295 mode = -1 or +1 (1/idim scaling on inverse [+1] if CFO_sclinv
296 was set in the call to csfft_scale_inverse_OMP)
297 idim = dimension (power of 2, possibly with a factor
298 of 3^p * 5^q also, for p,q=0..RMAX)
299 xc = input/output array
300 Automagically re-initializes itself when idim changes from
301 previous call. By AJ (Andrzej Jesmanowicz), modified by RWCox.
302 ----------------------------------------------------------------------*/
303
304 /*- Macro to do 1/N scaling on inverse:
305 if it is ordered by the user, and
306 if mode is positive, and
307 if this is not a recursive call. -*/
308
309 #define SCLINV \
310 if( CFO_sclinv[ithr] && mode > 0 && rec == 0 ){ \
311 register int qq ; register float ff = 1.0 / (float) idim ; \
312 for( qq=0 ; qq < idim ; qq++ ){ xc[qq].r *= ff ; xc[qq].i *= ff ; } }
313
csfft_cox_OMP(int mode,int idim,complex * xc)314 void csfft_cox_OMP( int mode , int idim , complex *xc )
315 {
316 DECLARE_ithr ;
317 register unsigned int m, n, i0, i1, i2, i3, k;
318 register complex *r0, *r1, *csp;
319 register float co, si, f0, f1, f2, f3, f4;
320 int rec = CFO_recC[ithr] ;
321
322 if( idim <= 1 ) return ; /* stoopid inpoot */
323
324 /*-- November 1998: maybe use the unrolled FFT routines --*/
325
326 /* INFO_message("csfft_cox_OMP(%d)",idim) ; */
327
328 switch( idim ){ /* none of these */
329 case 1: return; /* routines will */
330 case 2: fft2 (mode,xc); SCLINV; return; /* call csfft_cox */
331 case 4: fft4 (mode,xc); SCLINV; return; /* so don't need */
332 case 8: fft8 (mode,xc); SCLINV; return; /* to do rec++/-- */
333 case 16: fft16 (mode,xc); SCLINV; return;
334 case 32: fft32 (mode,xc); SCLINV; return;
335 case 64: fft64 (mode,xc); SCLINV; return;
336 #ifndef DONT_USE_4DEC
337 case 128: fft128 (mode,xc); SCLINV; return;
338 case 256: fft256 (mode,xc); SCLINV; return;
339 case 512: fft512 (mode,xc); SCLINV; return;
340 case 1024: fft1024 (mode,xc); SCLINV; return;
341 #if 1
342 case 2048: fft2048 (mode,xc); SCLINV; return;
343 case 4096: fft4096 (mode,xc); SCLINV; return;
344 case 8192: fft8192 (mode,xc); SCLINV; return;
345 case 16384: fft16384(mode,xc); SCLINV; return;
346 case 32768: fft32768(mode,xc); SCLINV; return;
347 case 65536: fft65536(mode,xc); SCLINV; return;
348 #endif
349 #endif
350 }
351
352 /* fall thru to here ==> not one of the shorter 2^p cases above! */
353
354 /*-- Aug 1999: deal with factors of 3 or 5 [might call csfft_cox_OMP again] --*/
355
356 if( idim%3 == 0 ){CFO_recC[ithr]++; fft_3dec_OMP(mode,idim,xc); CFO_recC[ithr]--; SCLINV; return;}
357 if( idim%5 == 0 ){CFO_recC[ithr]++; fft_5dec_OMP(mode,idim,xc); CFO_recC[ithr]--; SCLINV; return;}
358
359 /*-- below here: the old general power of 2 routine (i.e., for idim > 65536) --*/
360
361 /**-- perhaps initialize trig consts (for this thread) --**/
362
363 if( CFO_nold[ithr] != idim ) csfft_trigconsts( idim ) ;
364
365 /** Main loop starts here **/
366
367 n = idim;
368 i2 = idim >> 1;
369 i1 = 0;
370 csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */
371
372 /*-- swap data --*/
373
374 for (i0=0; i0 < n; i0 ++) {
375 if ( i1 > i0 ) {
376 r0 = xc + i0; r1 = xc + i1;
377 f1 = r0->r; f2 = r0->i;
378 r0->r = r1->r; r0->i = r1->i;
379 r1->r = f1; r1->i = f2;
380 }
381 m = i2;
382 while ( m && !(i1 < m) ) { i1 -= m; m >>= 1; }
383 i1 += m;
384 }
385
386 /*-- compute compute compute --*/
387
388 m = 1; k = 0;
389 while (n > m) {
390 i3 = m << 1;
391 for (i0=0; i0 < m; i0 ++) {
392 co = (csp + k)->r; si = (csp + k)->i;
393 for (i1=i0; i1 < n; i1 += i3) {
394 r0 = xc + i1; r1 = r0 + m;
395 f1 = r1->r * co; f2 = r1->i * si;
396 f3 = r1->r * si; f4 = r1->i * co;
397 f1 -= f2; f3 += f4;
398 f2 = r0->r; f4 = r0->i;
399 r1->r = f2 - f1; r1->i = f4 - f3;
400 r0->r = f2 + f1; r0->i = f4 + f3;
401 }
402 k++;
403 }
404 m = i3;
405 }
406
407 SCLINV ; return ;
408 }
409
410 /******************************************************************************
411 Routines below for unrolled FFT routines; generated with program fftprint.c.
412 *******************************************************************************/
413
414 /**************************************/
415 /* FFT routine unrolled of length 8 */
416
fft8(int mode,complex * xc)417 static void fft8( int mode , complex *xc )
418 {
419 DECLARE_ithr ;
420 register complex *csp , *xcx=xc;
421 register float f1,f2,f3,f4 ;
422
423 /* ININFO_message("fft8") ; */
424
425 /** perhaps initialize **/
426
427 if( CFO_nold[ithr] != 8 ) csfft_trigconsts( 8 ) ;
428
429 csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */
430
431 /** data swapping part **/
432
433 f1 = xcx[1].r ; f2 = xcx[1].i ;
434 xcx[1].r = xcx[4].r ; xcx[1].i = xcx[4].i ;
435 xcx[4].r = f1 ; xcx[4].i = f2 ;
436
437 f1 = xcx[3].r ; f2 = xcx[3].i ;
438 xcx[3].r = xcx[6].r ; xcx[3].i = xcx[6].i ;
439 xcx[6].r = f1 ; xcx[6].i = f2 ;
440
441 /** butterflying part **/
442
443 f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */
444 f2 = xcx[0].r ; f4 = xcx[0].i ;
445 xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ;
446 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
447
448 f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */
449 f2 = xcx[2].r ; f4 = xcx[2].i ;
450 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
451 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
452
453 f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */
454 f2 = xcx[4].r ; f4 = xcx[4].i ;
455 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
456 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
457
458 f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */
459 f2 = xcx[6].r ; f4 = xcx[6].i ;
460 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
461 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
462
463 f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */
464 f2 = xcx[0].r ; f4 = xcx[0].i ;
465 xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ;
466 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
467
468 f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */
469 f2 = xcx[4].r ; f4 = xcx[4].i ;
470 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
471 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
472
473 f1 = xcx[3].r * csp[2].r - xcx[3].i * csp[2].i ; /* twiddles */
474 f3 = xcx[3].r * csp[2].i + xcx[3].i * csp[2].r ;
475 f2 = xcx[1].r ; f4 = xcx[1].i ;
476 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
477 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
478
479 f1 = xcx[7].r * csp[2].r - xcx[7].i * csp[2].i ; /* twiddles */
480 f3 = xcx[7].r * csp[2].i + xcx[7].i * csp[2].r ;
481 f2 = xcx[5].r ; f4 = xcx[5].i ;
482 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
483 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
484
485 f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */
486 f2 = xcx[0].r ; f4 = xcx[0].i ;
487 xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ;
488 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
489
490 f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */
491 f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ;
492 f2 = xcx[1].r ; f4 = xcx[1].i ;
493 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
494 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
495
496 f1 = xcx[6].r * csp[5].r - xcx[6].i * csp[5].i ; /* twiddles */
497 f3 = xcx[6].r * csp[5].i + xcx[6].i * csp[5].r ;
498 f2 = xcx[2].r ; f4 = xcx[2].i ;
499 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
500 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
501
502 f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */
503 f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ;
504 f2 = xcx[3].r ; f4 = xcx[3].i ;
505 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
506 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
507
508 return ;
509 }
510
511 /**************************************/
512 /* FFT routine unrolled of length 16 */
513
fft16(int mode,complex * xc)514 static void fft16( int mode , complex *xc )
515 {
516 DECLARE_ithr ;
517 register complex *csp , *xcx=xc;
518 register float f1,f2,f3,f4 ;
519
520 /* ININFO_message("fft16") ; */
521
522 /** perhaps initialize **/
523
524 if( CFO_nold[ithr] != 16 ) csfft_trigconsts( 16 ) ;
525
526 csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */
527
528 /** data swapping part **/
529
530 f1 = xcx[1].r ; f2 = xcx[1].i ;
531 xcx[1].r = xcx[8].r ; xcx[1].i = xcx[8].i ;
532 xcx[8].r = f1 ; xcx[8].i = f2 ;
533
534 f1 = xcx[2].r ; f2 = xcx[2].i ;
535 xcx[2].r = xcx[4].r ; xcx[2].i = xcx[4].i ;
536 xcx[4].r = f1 ; xcx[4].i = f2 ;
537
538 f1 = xcx[3].r ; f2 = xcx[3].i ;
539 xcx[3].r = xcx[12].r ; xcx[3].i = xcx[12].i ;
540 xcx[12].r = f1 ; xcx[12].i = f2 ;
541
542 f1 = xcx[5].r ; f2 = xcx[5].i ;
543 xcx[5].r = xcx[10].r ; xcx[5].i = xcx[10].i ;
544 xcx[10].r = f1 ; xcx[10].i = f2 ;
545
546 f1 = xcx[7].r ; f2 = xcx[7].i ;
547 xcx[7].r = xcx[14].r ; xcx[7].i = xcx[14].i ;
548 xcx[14].r = f1 ; xcx[14].i = f2 ;
549
550 f1 = xcx[11].r ; f2 = xcx[11].i ;
551 xcx[11].r = xcx[13].r ; xcx[11].i = xcx[13].i ;
552 xcx[13].r = f1 ; xcx[13].i = f2 ;
553
554 /** butterflying part **/
555
556 f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */
557 f2 = xcx[0].r ; f4 = xcx[0].i ;
558 xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ;
559 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
560
561 f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */
562 f2 = xcx[2].r ; f4 = xcx[2].i ;
563 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
564 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
565
566 f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */
567 f2 = xcx[4].r ; f4 = xcx[4].i ;
568 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
569 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
570
571 f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */
572 f2 = xcx[6].r ; f4 = xcx[6].i ;
573 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
574 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
575
576 f1 = xcx[9].r ; f3 = xcx[9].i ; /* cos=1 sin=0 */
577 f2 = xcx[8].r ; f4 = xcx[8].i ;
578 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
579 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
580
581 f1 = xcx[11].r ; f3 = xcx[11].i ; /* cos=1 sin=0 */
582 f2 = xcx[10].r ; f4 = xcx[10].i ;
583 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
584 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
585
586 f1 = xcx[13].r ; f3 = xcx[13].i ; /* cos=1 sin=0 */
587 f2 = xcx[12].r ; f4 = xcx[12].i ;
588 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
589 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
590
591 f1 = xcx[15].r ; f3 = xcx[15].i ; /* cos=1 sin=0 */
592 f2 = xcx[14].r ; f4 = xcx[14].i ;
593 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
594 xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ;
595
596 f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */
597 f2 = xcx[0].r ; f4 = xcx[0].i ;
598 xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ;
599 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
600
601 f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */
602 f2 = xcx[4].r ; f4 = xcx[4].i ;
603 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
604 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
605
606 f1 = xcx[10].r ; f3 = xcx[10].i ; /* cos=1 sin=0 */
607 f2 = xcx[8].r ; f4 = xcx[8].i ;
608 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
609 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
610
611 f1 = xcx[14].r ; f3 = xcx[14].i ; /* cos=1 sin=0 */
612 f2 = xcx[12].r ; f4 = xcx[12].i ;
613 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
614 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
615
616 f1 = xcx[3].r * csp[2].r - xcx[3].i * csp[2].i ; /* twiddles */
617 f3 = xcx[3].r * csp[2].i + xcx[3].i * csp[2].r ;
618 f2 = xcx[1].r ; f4 = xcx[1].i ;
619 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
620 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
621
622 f1 = xcx[7].r * csp[2].r - xcx[7].i * csp[2].i ; /* twiddles */
623 f3 = xcx[7].r * csp[2].i + xcx[7].i * csp[2].r ;
624 f2 = xcx[5].r ; f4 = xcx[5].i ;
625 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
626 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
627
628 f1 = xcx[11].r * csp[2].r - xcx[11].i * csp[2].i ; /* twiddles */
629 f3 = xcx[11].r * csp[2].i + xcx[11].i * csp[2].r ;
630 f2 = xcx[9].r ; f4 = xcx[9].i ;
631 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
632 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
633
634 f1 = xcx[15].r * csp[2].r - xcx[15].i * csp[2].i ; /* twiddles */
635 f3 = xcx[15].r * csp[2].i + xcx[15].i * csp[2].r ;
636 f2 = xcx[13].r ; f4 = xcx[13].i ;
637 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
638 xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ;
639
640 f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */
641 f2 = xcx[0].r ; f4 = xcx[0].i ;
642 xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ;
643 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
644
645 f1 = xcx[12].r ; f3 = xcx[12].i ; /* cos=1 sin=0 */
646 f2 = xcx[8].r ; f4 = xcx[8].i ;
647 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
648 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
649
650 f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */
651 f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ;
652 f2 = xcx[1].r ; f4 = xcx[1].i ;
653 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
654 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
655
656 f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */
657 f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ;
658 f2 = xcx[9].r ; f4 = xcx[9].i ;
659 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
660 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
661
662 f1 = xcx[6].r * csp[5].r - xcx[6].i * csp[5].i ; /* twiddles */
663 f3 = xcx[6].r * csp[5].i + xcx[6].i * csp[5].r ;
664 f2 = xcx[2].r ; f4 = xcx[2].i ;
665 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
666 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
667
668 f1 = xcx[14].r * csp[5].r - xcx[14].i * csp[5].i ; /* twiddles */
669 f3 = xcx[14].r * csp[5].i + xcx[14].i * csp[5].r ;
670 f2 = xcx[10].r ; f4 = xcx[10].i ;
671 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
672 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
673
674 f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */
675 f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ;
676 f2 = xcx[3].r ; f4 = xcx[3].i ;
677 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
678 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
679
680 f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */
681 f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ;
682 f2 = xcx[11].r ; f4 = xcx[11].i ;
683 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
684 xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ;
685
686 f1 = xcx[8].r ; f3 = xcx[8].i ; /* cos=1 sin=0 */
687 f2 = xcx[0].r ; f4 = xcx[0].i ;
688 xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ;
689 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
690
691 f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */
692 f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ;
693 f2 = xcx[1].r ; f4 = xcx[1].i ;
694 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
695 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
696
697 f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */
698 f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ;
699 f2 = xcx[2].r ; f4 = xcx[2].i ;
700 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
701 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
702
703 f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */
704 f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ;
705 f2 = xcx[3].r ; f4 = xcx[3].i ;
706 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
707 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
708
709 f1 = xcx[12].r * csp[11].r - xcx[12].i * csp[11].i ; /* twiddles */
710 f3 = xcx[12].r * csp[11].i + xcx[12].i * csp[11].r ;
711 f2 = xcx[4].r ; f4 = xcx[4].i ;
712 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
713 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
714
715 f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */
716 f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ;
717 f2 = xcx[5].r ; f4 = xcx[5].i ;
718 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
719 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
720
721 f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */
722 f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ;
723 f2 = xcx[6].r ; f4 = xcx[6].i ;
724 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
725 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
726
727 f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */
728 f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ;
729 f2 = xcx[7].r ; f4 = xcx[7].i ;
730 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
731 xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ;
732
733 return ;
734 }
735
736 /**************************************/
737 /* FFT routine unrolled of length 32 */
738
fft32(int mode,complex * xc)739 static void fft32( int mode , complex *xc )
740 {
741 DECLARE_ithr ;
742 register complex *csp , *xcx=xc;
743 register float f1,f2,f3,f4 ;
744
745 /* ININFO_message("fft32") ; */
746
747 /** perhaps initialize **/
748
749 if( CFO_nold[ithr] != 32 ) csfft_trigconsts( 32 ) ;
750
751 csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */
752
753 /** data swapping part **/
754
755 f1 = xcx[1].r ; f2 = xcx[1].i ;
756 xcx[1].r = xcx[16].r ; xcx[1].i = xcx[16].i ;
757 xcx[16].r = f1 ; xcx[16].i = f2 ;
758
759 f1 = xcx[2].r ; f2 = xcx[2].i ;
760 xcx[2].r = xcx[8].r ; xcx[2].i = xcx[8].i ;
761 xcx[8].r = f1 ; xcx[8].i = f2 ;
762
763 f1 = xcx[3].r ; f2 = xcx[3].i ;
764 xcx[3].r = xcx[24].r ; xcx[3].i = xcx[24].i ;
765 xcx[24].r = f1 ; xcx[24].i = f2 ;
766
767 f1 = xcx[5].r ; f2 = xcx[5].i ;
768 xcx[5].r = xcx[20].r ; xcx[5].i = xcx[20].i ;
769 xcx[20].r = f1 ; xcx[20].i = f2 ;
770
771 f1 = xcx[6].r ; f2 = xcx[6].i ;
772 xcx[6].r = xcx[12].r ; xcx[6].i = xcx[12].i ;
773 xcx[12].r = f1 ; xcx[12].i = f2 ;
774
775 f1 = xcx[7].r ; f2 = xcx[7].i ;
776 xcx[7].r = xcx[28].r ; xcx[7].i = xcx[28].i ;
777 xcx[28].r = f1 ; xcx[28].i = f2 ;
778
779 f1 = xcx[9].r ; f2 = xcx[9].i ;
780 xcx[9].r = xcx[18].r ; xcx[9].i = xcx[18].i ;
781 xcx[18].r = f1 ; xcx[18].i = f2 ;
782
783 f1 = xcx[11].r ; f2 = xcx[11].i ;
784 xcx[11].r = xcx[26].r ; xcx[11].i = xcx[26].i ;
785 xcx[26].r = f1 ; xcx[26].i = f2 ;
786
787 f1 = xcx[13].r ; f2 = xcx[13].i ;
788 xcx[13].r = xcx[22].r ; xcx[13].i = xcx[22].i ;
789 xcx[22].r = f1 ; xcx[22].i = f2 ;
790
791 f1 = xcx[15].r ; f2 = xcx[15].i ;
792 xcx[15].r = xcx[30].r ; xcx[15].i = xcx[30].i ;
793 xcx[30].r = f1 ; xcx[30].i = f2 ;
794
795 f1 = xcx[19].r ; f2 = xcx[19].i ;
796 xcx[19].r = xcx[25].r ; xcx[19].i = xcx[25].i ;
797 xcx[25].r = f1 ; xcx[25].i = f2 ;
798
799 f1 = xcx[23].r ; f2 = xcx[23].i ;
800 xcx[23].r = xcx[29].r ; xcx[23].i = xcx[29].i ;
801 xcx[29].r = f1 ; xcx[29].i = f2 ;
802
803 /** butterflying part **/
804
805 f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */
806 f2 = xcx[0].r ; f4 = xcx[0].i ;
807 xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ;
808 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
809
810 f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */
811 f2 = xcx[2].r ; f4 = xcx[2].i ;
812 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
813 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
814
815 f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */
816 f2 = xcx[4].r ; f4 = xcx[4].i ;
817 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
818 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
819
820 f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */
821 f2 = xcx[6].r ; f4 = xcx[6].i ;
822 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
823 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
824
825 f1 = xcx[9].r ; f3 = xcx[9].i ; /* cos=1 sin=0 */
826 f2 = xcx[8].r ; f4 = xcx[8].i ;
827 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
828 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
829
830 f1 = xcx[11].r ; f3 = xcx[11].i ; /* cos=1 sin=0 */
831 f2 = xcx[10].r ; f4 = xcx[10].i ;
832 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
833 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
834
835 f1 = xcx[13].r ; f3 = xcx[13].i ; /* cos=1 sin=0 */
836 f2 = xcx[12].r ; f4 = xcx[12].i ;
837 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
838 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
839
840 f1 = xcx[15].r ; f3 = xcx[15].i ; /* cos=1 sin=0 */
841 f2 = xcx[14].r ; f4 = xcx[14].i ;
842 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
843 xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ;
844
845 f1 = xcx[17].r ; f3 = xcx[17].i ; /* cos=1 sin=0 */
846 f2 = xcx[16].r ; f4 = xcx[16].i ;
847 xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ;
848 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
849
850 f1 = xcx[19].r ; f3 = xcx[19].i ; /* cos=1 sin=0 */
851 f2 = xcx[18].r ; f4 = xcx[18].i ;
852 xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ;
853 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
854
855 f1 = xcx[21].r ; f3 = xcx[21].i ; /* cos=1 sin=0 */
856 f2 = xcx[20].r ; f4 = xcx[20].i ;
857 xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ;
858 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
859
860 f1 = xcx[23].r ; f3 = xcx[23].i ; /* cos=1 sin=0 */
861 f2 = xcx[22].r ; f4 = xcx[22].i ;
862 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
863 xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ;
864
865 f1 = xcx[25].r ; f3 = xcx[25].i ; /* cos=1 sin=0 */
866 f2 = xcx[24].r ; f4 = xcx[24].i ;
867 xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ;
868 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
869
870 f1 = xcx[27].r ; f3 = xcx[27].i ; /* cos=1 sin=0 */
871 f2 = xcx[26].r ; f4 = xcx[26].i ;
872 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
873 xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ;
874
875 f1 = xcx[29].r ; f3 = xcx[29].i ; /* cos=1 sin=0 */
876 f2 = xcx[28].r ; f4 = xcx[28].i ;
877 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
878 xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ;
879
880 f1 = xcx[31].r ; f3 = xcx[31].i ; /* cos=1 sin=0 */
881 f2 = xcx[30].r ; f4 = xcx[30].i ;
882 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
883 xcx[30].r = f2+f1 ; xcx[30].i = f4+f3 ;
884
885 f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */
886 f2 = xcx[0].r ; f4 = xcx[0].i ;
887 xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ;
888 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
889
890 f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */
891 f2 = xcx[4].r ; f4 = xcx[4].i ;
892 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
893 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
894
895 f1 = xcx[10].r ; f3 = xcx[10].i ; /* cos=1 sin=0 */
896 f2 = xcx[8].r ; f4 = xcx[8].i ;
897 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
898 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
899
900 f1 = xcx[14].r ; f3 = xcx[14].i ; /* cos=1 sin=0 */
901 f2 = xcx[12].r ; f4 = xcx[12].i ;
902 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
903 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
904
905 f1 = xcx[18].r ; f3 = xcx[18].i ; /* cos=1 sin=0 */
906 f2 = xcx[16].r ; f4 = xcx[16].i ;
907 xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ;
908 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
909
910 f1 = xcx[22].r ; f3 = xcx[22].i ; /* cos=1 sin=0 */
911 f2 = xcx[20].r ; f4 = xcx[20].i ;
912 xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ;
913 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
914
915 f1 = xcx[26].r ; f3 = xcx[26].i ; /* cos=1 sin=0 */
916 f2 = xcx[24].r ; f4 = xcx[24].i ;
917 xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ;
918 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
919
920 f1 = xcx[30].r ; f3 = xcx[30].i ; /* cos=1 sin=0 */
921 f2 = xcx[28].r ; f4 = xcx[28].i ;
922 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
923 xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ;
924
925 f1 = - xcx[3].i * csp[2].i ; /* cos=0 twiddles */
926 f3 = xcx[3].r * csp[2].i ;
927 f2 = xcx[1].r ; f4 = xcx[1].i ;
928 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
929 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
930
931 f1 = - xcx[7].i * csp[2].i ; /* cos=0 twiddles */
932 f3 = xcx[7].r * csp[2].i ;
933 f2 = xcx[5].r ; f4 = xcx[5].i ;
934 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
935 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
936
937 f1 = - xcx[11].i * csp[2].i ; /* cos=0 twiddles */
938 f3 = xcx[11].r * csp[2].i ;
939 f2 = xcx[9].r ; f4 = xcx[9].i ;
940 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
941 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
942
943 f1 = - xcx[15].i * csp[2].i ; /* cos=0 twiddles */
944 f3 = xcx[15].r * csp[2].i ;
945 f2 = xcx[13].r ; f4 = xcx[13].i ;
946 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
947 xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ;
948
949 f1 = - xcx[19].i * csp[2].i ; /* cos=0 twiddles */
950 f3 = xcx[19].r * csp[2].i ;
951 f2 = xcx[17].r ; f4 = xcx[17].i ;
952 xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ;
953 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
954
955 f1 = - xcx[23].i * csp[2].i ; /* cos=0 twiddles */
956 f3 = xcx[23].r * csp[2].i ;
957 f2 = xcx[21].r ; f4 = xcx[21].i ;
958 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
959 xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ;
960
961 f1 = - xcx[27].i * csp[2].i ; /* cos=0 twiddles */
962 f3 = xcx[27].r * csp[2].i ;
963 f2 = xcx[25].r ; f4 = xcx[25].i ;
964 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
965 xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ;
966
967 f1 = - xcx[31].i * csp[2].i ; /* cos=0 twiddles */
968 f3 = xcx[31].r * csp[2].i ;
969 f2 = xcx[29].r ; f4 = xcx[29].i ;
970 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
971 xcx[29].r = f2+f1 ; xcx[29].i = f4+f3 ;
972
973 f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */
974 f2 = xcx[0].r ; f4 = xcx[0].i ;
975 xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ;
976 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
977
978 f1 = xcx[12].r ; f3 = xcx[12].i ; /* cos=1 sin=0 */
979 f2 = xcx[8].r ; f4 = xcx[8].i ;
980 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
981 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
982
983 f1 = xcx[20].r ; f3 = xcx[20].i ; /* cos=1 sin=0 */
984 f2 = xcx[16].r ; f4 = xcx[16].i ;
985 xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ;
986 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
987
988 f1 = xcx[28].r ; f3 = xcx[28].i ; /* cos=1 sin=0 */
989 f2 = xcx[24].r ; f4 = xcx[24].i ;
990 xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ;
991 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
992
993 f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */
994 f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ;
995 f2 = xcx[1].r ; f4 = xcx[1].i ;
996 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
997 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
998
999 f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */
1000 f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ;
1001 f2 = xcx[9].r ; f4 = xcx[9].i ;
1002 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
1003 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
1004
1005 f1 = xcx[21].r * csp[4].r - xcx[21].i * csp[4].i ; /* twiddles */
1006 f3 = xcx[21].r * csp[4].i + xcx[21].i * csp[4].r ;
1007 f2 = xcx[17].r ; f4 = xcx[17].i ;
1008 xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ;
1009 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
1010
1011 f1 = xcx[29].r * csp[4].r - xcx[29].i * csp[4].i ; /* twiddles */
1012 f3 = xcx[29].r * csp[4].i + xcx[29].i * csp[4].r ;
1013 f2 = xcx[25].r ; f4 = xcx[25].i ;
1014 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
1015 xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ;
1016
1017 f1 = - xcx[6].i * csp[5].i ; /* cos=0 twiddles */
1018 f3 = xcx[6].r * csp[5].i ;
1019 f2 = xcx[2].r ; f4 = xcx[2].i ;
1020 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
1021 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
1022
1023 f1 = - xcx[14].i * csp[5].i ; /* cos=0 twiddles */
1024 f3 = xcx[14].r * csp[5].i ;
1025 f2 = xcx[10].r ; f4 = xcx[10].i ;
1026 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
1027 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
1028
1029 f1 = - xcx[22].i * csp[5].i ; /* cos=0 twiddles */
1030 f3 = xcx[22].r * csp[5].i ;
1031 f2 = xcx[18].r ; f4 = xcx[18].i ;
1032 xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ;
1033 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
1034
1035 f1 = - xcx[30].i * csp[5].i ; /* cos=0 twiddles */
1036 f3 = xcx[30].r * csp[5].i ;
1037 f2 = xcx[26].r ; f4 = xcx[26].i ;
1038 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
1039 xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ;
1040
1041 f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */
1042 f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ;
1043 f2 = xcx[3].r ; f4 = xcx[3].i ;
1044 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
1045 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
1046
1047 f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */
1048 f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ;
1049 f2 = xcx[11].r ; f4 = xcx[11].i ;
1050 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
1051 xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ;
1052
1053 f1 = xcx[23].r * csp[6].r - xcx[23].i * csp[6].i ; /* twiddles */
1054 f3 = xcx[23].r * csp[6].i + xcx[23].i * csp[6].r ;
1055 f2 = xcx[19].r ; f4 = xcx[19].i ;
1056 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
1057 xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ;
1058
1059 f1 = xcx[31].r * csp[6].r - xcx[31].i * csp[6].i ; /* twiddles */
1060 f3 = xcx[31].r * csp[6].i + xcx[31].i * csp[6].r ;
1061 f2 = xcx[27].r ; f4 = xcx[27].i ;
1062 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
1063 xcx[27].r = f2+f1 ; xcx[27].i = f4+f3 ;
1064
1065 f1 = xcx[8].r ; f3 = xcx[8].i ; /* cos=1 sin=0 */
1066 f2 = xcx[0].r ; f4 = xcx[0].i ;
1067 xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ;
1068 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
1069
1070 f1 = xcx[24].r ; f3 = xcx[24].i ; /* cos=1 sin=0 */
1071 f2 = xcx[16].r ; f4 = xcx[16].i ;
1072 xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ;
1073 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
1074
1075 f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */
1076 f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ;
1077 f2 = xcx[1].r ; f4 = xcx[1].i ;
1078 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
1079 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
1080
1081 f1 = xcx[25].r * csp[8].r - xcx[25].i * csp[8].i ; /* twiddles */
1082 f3 = xcx[25].r * csp[8].i + xcx[25].i * csp[8].r ;
1083 f2 = xcx[17].r ; f4 = xcx[17].i ;
1084 xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ;
1085 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
1086
1087 f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */
1088 f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ;
1089 f2 = xcx[2].r ; f4 = xcx[2].i ;
1090 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
1091 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
1092
1093 f1 = xcx[26].r * csp[9].r - xcx[26].i * csp[9].i ; /* twiddles */
1094 f3 = xcx[26].r * csp[9].i + xcx[26].i * csp[9].r ;
1095 f2 = xcx[18].r ; f4 = xcx[18].i ;
1096 xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ;
1097 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
1098
1099 f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */
1100 f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ;
1101 f2 = xcx[3].r ; f4 = xcx[3].i ;
1102 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
1103 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
1104
1105 f1 = xcx[27].r * csp[10].r - xcx[27].i * csp[10].i ; /* twiddles */
1106 f3 = xcx[27].r * csp[10].i + xcx[27].i * csp[10].r ;
1107 f2 = xcx[19].r ; f4 = xcx[19].i ;
1108 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
1109 xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ;
1110
1111 f1 = - xcx[12].i * csp[11].i ; /* cos=0 twiddles */
1112 f3 = xcx[12].r * csp[11].i ;
1113 f2 = xcx[4].r ; f4 = xcx[4].i ;
1114 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
1115 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
1116
1117 f1 = - xcx[28].i * csp[11].i ; /* cos=0 twiddles */
1118 f3 = xcx[28].r * csp[11].i ;
1119 f2 = xcx[20].r ; f4 = xcx[20].i ;
1120 xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ;
1121 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
1122
1123 f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */
1124 f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ;
1125 f2 = xcx[5].r ; f4 = xcx[5].i ;
1126 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
1127 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
1128
1129 f1 = xcx[29].r * csp[12].r - xcx[29].i * csp[12].i ; /* twiddles */
1130 f3 = xcx[29].r * csp[12].i + xcx[29].i * csp[12].r ;
1131 f2 = xcx[21].r ; f4 = xcx[21].i ;
1132 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
1133 xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ;
1134
1135 f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */
1136 f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ;
1137 f2 = xcx[6].r ; f4 = xcx[6].i ;
1138 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
1139 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
1140
1141 f1 = xcx[30].r * csp[13].r - xcx[30].i * csp[13].i ; /* twiddles */
1142 f3 = xcx[30].r * csp[13].i + xcx[30].i * csp[13].r ;
1143 f2 = xcx[22].r ; f4 = xcx[22].i ;
1144 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
1145 xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ;
1146
1147 f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */
1148 f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ;
1149 f2 = xcx[7].r ; f4 = xcx[7].i ;
1150 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
1151 xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ;
1152
1153 f1 = xcx[31].r * csp[14].r - xcx[31].i * csp[14].i ; /* twiddles */
1154 f3 = xcx[31].r * csp[14].i + xcx[31].i * csp[14].r ;
1155 f2 = xcx[23].r ; f4 = xcx[23].i ;
1156 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
1157 xcx[23].r = f2+f1 ; xcx[23].i = f4+f3 ;
1158
1159 f1 = xcx[16].r ; f3 = xcx[16].i ; /* cos=1 sin=0 */
1160 f2 = xcx[0].r ; f4 = xcx[0].i ;
1161 xcx[16].r = f2-f1 ; xcx[16].i = f4-f3 ;
1162 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
1163
1164 f1 = xcx[17].r * csp[16].r - xcx[17].i * csp[16].i ; /* twiddles */
1165 f3 = xcx[17].r * csp[16].i + xcx[17].i * csp[16].r ;
1166 f2 = xcx[1].r ; f4 = xcx[1].i ;
1167 xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ;
1168 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
1169
1170 f1 = xcx[18].r * csp[17].r - xcx[18].i * csp[17].i ; /* twiddles */
1171 f3 = xcx[18].r * csp[17].i + xcx[18].i * csp[17].r ;
1172 f2 = xcx[2].r ; f4 = xcx[2].i ;
1173 xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ;
1174 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
1175
1176 f1 = xcx[19].r * csp[18].r - xcx[19].i * csp[18].i ; /* twiddles */
1177 f3 = xcx[19].r * csp[18].i + xcx[19].i * csp[18].r ;
1178 f2 = xcx[3].r ; f4 = xcx[3].i ;
1179 xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ;
1180 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
1181
1182 f1 = xcx[20].r * csp[19].r - xcx[20].i * csp[19].i ; /* twiddles */
1183 f3 = xcx[20].r * csp[19].i + xcx[20].i * csp[19].r ;
1184 f2 = xcx[4].r ; f4 = xcx[4].i ;
1185 xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ;
1186 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
1187
1188 f1 = xcx[21].r * csp[20].r - xcx[21].i * csp[20].i ; /* twiddles */
1189 f3 = xcx[21].r * csp[20].i + xcx[21].i * csp[20].r ;
1190 f2 = xcx[5].r ; f4 = xcx[5].i ;
1191 xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ;
1192 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
1193
1194 f1 = xcx[22].r * csp[21].r - xcx[22].i * csp[21].i ; /* twiddles */
1195 f3 = xcx[22].r * csp[21].i + xcx[22].i * csp[21].r ;
1196 f2 = xcx[6].r ; f4 = xcx[6].i ;
1197 xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ;
1198 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
1199
1200 f1 = xcx[23].r * csp[22].r - xcx[23].i * csp[22].i ; /* twiddles */
1201 f3 = xcx[23].r * csp[22].i + xcx[23].i * csp[22].r ;
1202 f2 = xcx[7].r ; f4 = xcx[7].i ;
1203 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
1204 xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ;
1205
1206 f1 = - xcx[24].i * csp[23].i ; /* cos=0 twiddles */
1207 f3 = xcx[24].r * csp[23].i ;
1208 f2 = xcx[8].r ; f4 = xcx[8].i ;
1209 xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ;
1210 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
1211
1212 f1 = xcx[25].r * csp[24].r - xcx[25].i * csp[24].i ; /* twiddles */
1213 f3 = xcx[25].r * csp[24].i + xcx[25].i * csp[24].r ;
1214 f2 = xcx[9].r ; f4 = xcx[9].i ;
1215 xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ;
1216 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
1217
1218 f1 = xcx[26].r * csp[25].r - xcx[26].i * csp[25].i ; /* twiddles */
1219 f3 = xcx[26].r * csp[25].i + xcx[26].i * csp[25].r ;
1220 f2 = xcx[10].r ; f4 = xcx[10].i ;
1221 xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ;
1222 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
1223
1224 f1 = xcx[27].r * csp[26].r - xcx[27].i * csp[26].i ; /* twiddles */
1225 f3 = xcx[27].r * csp[26].i + xcx[27].i * csp[26].r ;
1226 f2 = xcx[11].r ; f4 = xcx[11].i ;
1227 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
1228 xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ;
1229
1230 f1 = xcx[28].r * csp[27].r - xcx[28].i * csp[27].i ; /* twiddles */
1231 f3 = xcx[28].r * csp[27].i + xcx[28].i * csp[27].r ;
1232 f2 = xcx[12].r ; f4 = xcx[12].i ;
1233 xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ;
1234 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
1235
1236 f1 = xcx[29].r * csp[28].r - xcx[29].i * csp[28].i ; /* twiddles */
1237 f3 = xcx[29].r * csp[28].i + xcx[29].i * csp[28].r ;
1238 f2 = xcx[13].r ; f4 = xcx[13].i ;
1239 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
1240 xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ;
1241
1242 f1 = xcx[30].r * csp[29].r - xcx[30].i * csp[29].i ; /* twiddles */
1243 f3 = xcx[30].r * csp[29].i + xcx[30].i * csp[29].r ;
1244 f2 = xcx[14].r ; f4 = xcx[14].i ;
1245 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
1246 xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ;
1247
1248 f1 = xcx[31].r * csp[30].r - xcx[31].i * csp[30].i ; /* twiddles */
1249 f3 = xcx[31].r * csp[30].i + xcx[31].i * csp[30].r ;
1250 f2 = xcx[15].r ; f4 = xcx[15].i ;
1251 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
1252 xcx[15].r = f2+f1 ; xcx[15].i = f4+f3 ;
1253
1254 return ;
1255 }
1256
1257 /**************************************/
1258 /* FFT routine unrolled of length 64 */
1259
fft64(int mode,complex * xc)1260 void fft64( int mode , complex * xc )
1261 {
1262 DECLARE_ithr ;
1263 register complex *csp , *xcx=xc;
1264 register float f1,f2,f3,f4 ;
1265
1266 /* ININFO_message("fft64") ; */
1267
1268 /** perhaps initialize **/
1269
1270 if( CFO_nold[ithr] != 64 ) csfft_trigconsts( 64 ) ;
1271
1272 csp = (mode > 0) ? CFO_csplus[ithr] : CFO_csminus[ithr] ; /* choose const array */
1273
1274 /** data swapping part **/
1275
1276 f1 = xcx[1].r ; f2 = xcx[1].i ;
1277 xcx[1].r = xcx[32].r ; xcx[1].i = xcx[32].i ;
1278 xcx[32].r = f1 ; xcx[32].i = f2 ;
1279
1280 f1 = xcx[2].r ; f2 = xcx[2].i ;
1281 xcx[2].r = xcx[16].r ; xcx[2].i = xcx[16].i ;
1282 xcx[16].r = f1 ; xcx[16].i = f2 ;
1283
1284 f1 = xcx[3].r ; f2 = xcx[3].i ;
1285 xcx[3].r = xcx[48].r ; xcx[3].i = xcx[48].i ;
1286 xcx[48].r = f1 ; xcx[48].i = f2 ;
1287
1288 f1 = xcx[4].r ; f2 = xcx[4].i ;
1289 xcx[4].r = xcx[8].r ; xcx[4].i = xcx[8].i ;
1290 xcx[8].r = f1 ; xcx[8].i = f2 ;
1291
1292 f1 = xcx[5].r ; f2 = xcx[5].i ;
1293 xcx[5].r = xcx[40].r ; xcx[5].i = xcx[40].i ;
1294 xcx[40].r = f1 ; xcx[40].i = f2 ;
1295
1296 f1 = xcx[6].r ; f2 = xcx[6].i ;
1297 xcx[6].r = xcx[24].r ; xcx[6].i = xcx[24].i ;
1298 xcx[24].r = f1 ; xcx[24].i = f2 ;
1299
1300 f1 = xcx[7].r ; f2 = xcx[7].i ;
1301 xcx[7].r = xcx[56].r ; xcx[7].i = xcx[56].i ;
1302 xcx[56].r = f1 ; xcx[56].i = f2 ;
1303
1304 f1 = xcx[9].r ; f2 = xcx[9].i ;
1305 xcx[9].r = xcx[36].r ; xcx[9].i = xcx[36].i ;
1306 xcx[36].r = f1 ; xcx[36].i = f2 ;
1307
1308 f1 = xcx[10].r ; f2 = xcx[10].i ;
1309 xcx[10].r = xcx[20].r ; xcx[10].i = xcx[20].i ;
1310 xcx[20].r = f1 ; xcx[20].i = f2 ;
1311
1312 f1 = xcx[11].r ; f2 = xcx[11].i ;
1313 xcx[11].r = xcx[52].r ; xcx[11].i = xcx[52].i ;
1314 xcx[52].r = f1 ; xcx[52].i = f2 ;
1315
1316 f1 = xcx[13].r ; f2 = xcx[13].i ;
1317 xcx[13].r = xcx[44].r ; xcx[13].i = xcx[44].i ;
1318 xcx[44].r = f1 ; xcx[44].i = f2 ;
1319
1320 f1 = xcx[14].r ; f2 = xcx[14].i ;
1321 xcx[14].r = xcx[28].r ; xcx[14].i = xcx[28].i ;
1322 xcx[28].r = f1 ; xcx[28].i = f2 ;
1323
1324 f1 = xcx[15].r ; f2 = xcx[15].i ;
1325 xcx[15].r = xcx[60].r ; xcx[15].i = xcx[60].i ;
1326 xcx[60].r = f1 ; xcx[60].i = f2 ;
1327
1328 f1 = xcx[17].r ; f2 = xcx[17].i ;
1329 xcx[17].r = xcx[34].r ; xcx[17].i = xcx[34].i ;
1330 xcx[34].r = f1 ; xcx[34].i = f2 ;
1331
1332 f1 = xcx[19].r ; f2 = xcx[19].i ;
1333 xcx[19].r = xcx[50].r ; xcx[19].i = xcx[50].i ;
1334 xcx[50].r = f1 ; xcx[50].i = f2 ;
1335
1336 f1 = xcx[21].r ; f2 = xcx[21].i ;
1337 xcx[21].r = xcx[42].r ; xcx[21].i = xcx[42].i ;
1338 xcx[42].r = f1 ; xcx[42].i = f2 ;
1339
1340 f1 = xcx[22].r ; f2 = xcx[22].i ;
1341 xcx[22].r = xcx[26].r ; xcx[22].i = xcx[26].i ;
1342 xcx[26].r = f1 ; xcx[26].i = f2 ;
1343
1344 f1 = xcx[23].r ; f2 = xcx[23].i ;
1345 xcx[23].r = xcx[58].r ; xcx[23].i = xcx[58].i ;
1346 xcx[58].r = f1 ; xcx[58].i = f2 ;
1347
1348 f1 = xcx[25].r ; f2 = xcx[25].i ;
1349 xcx[25].r = xcx[38].r ; xcx[25].i = xcx[38].i ;
1350 xcx[38].r = f1 ; xcx[38].i = f2 ;
1351
1352 f1 = xcx[27].r ; f2 = xcx[27].i ;
1353 xcx[27].r = xcx[54].r ; xcx[27].i = xcx[54].i ;
1354 xcx[54].r = f1 ; xcx[54].i = f2 ;
1355
1356 f1 = xcx[29].r ; f2 = xcx[29].i ;
1357 xcx[29].r = xcx[46].r ; xcx[29].i = xcx[46].i ;
1358 xcx[46].r = f1 ; xcx[46].i = f2 ;
1359
1360 f1 = xcx[31].r ; f2 = xcx[31].i ;
1361 xcx[31].r = xcx[62].r ; xcx[31].i = xcx[62].i ;
1362 xcx[62].r = f1 ; xcx[62].i = f2 ;
1363
1364 f1 = xcx[35].r ; f2 = xcx[35].i ;
1365 xcx[35].r = xcx[49].r ; xcx[35].i = xcx[49].i ;
1366 xcx[49].r = f1 ; xcx[49].i = f2 ;
1367
1368 f1 = xcx[37].r ; f2 = xcx[37].i ;
1369 xcx[37].r = xcx[41].r ; xcx[37].i = xcx[41].i ;
1370 xcx[41].r = f1 ; xcx[41].i = f2 ;
1371
1372 f1 = xcx[39].r ; f2 = xcx[39].i ;
1373 xcx[39].r = xcx[57].r ; xcx[39].i = xcx[57].i ;
1374 xcx[57].r = f1 ; xcx[57].i = f2 ;
1375
1376 f1 = xcx[43].r ; f2 = xcx[43].i ;
1377 xcx[43].r = xcx[53].r ; xcx[43].i = xcx[53].i ;
1378 xcx[53].r = f1 ; xcx[53].i = f2 ;
1379
1380 f1 = xcx[47].r ; f2 = xcx[47].i ;
1381 xcx[47].r = xcx[61].r ; xcx[47].i = xcx[61].i ;
1382 xcx[61].r = f1 ; xcx[61].i = f2 ;
1383
1384 f1 = xcx[55].r ; f2 = xcx[55].i ;
1385 xcx[55].r = xcx[59].r ; xcx[55].i = xcx[59].i ;
1386 xcx[59].r = f1 ; xcx[59].i = f2 ;
1387
1388 /** butterflying part **/
1389
1390 f1 = xcx[1].r ; f3 = xcx[1].i ; /* cos=1 sin=0 */
1391 f2 = xcx[0].r ; f4 = xcx[0].i ;
1392 xcx[1].r = f2-f1 ; xcx[1].i = f4-f3 ;
1393 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
1394
1395 f1 = xcx[3].r ; f3 = xcx[3].i ; /* cos=1 sin=0 */
1396 f2 = xcx[2].r ; f4 = xcx[2].i ;
1397 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
1398 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
1399
1400 f1 = xcx[5].r ; f3 = xcx[5].i ; /* cos=1 sin=0 */
1401 f2 = xcx[4].r ; f4 = xcx[4].i ;
1402 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
1403 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
1404
1405 f1 = xcx[7].r ; f3 = xcx[7].i ; /* cos=1 sin=0 */
1406 f2 = xcx[6].r ; f4 = xcx[6].i ;
1407 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
1408 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
1409
1410 f1 = xcx[9].r ; f3 = xcx[9].i ; /* cos=1 sin=0 */
1411 f2 = xcx[8].r ; f4 = xcx[8].i ;
1412 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
1413 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
1414
1415 f1 = xcx[11].r ; f3 = xcx[11].i ; /* cos=1 sin=0 */
1416 f2 = xcx[10].r ; f4 = xcx[10].i ;
1417 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
1418 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
1419
1420 f1 = xcx[13].r ; f3 = xcx[13].i ; /* cos=1 sin=0 */
1421 f2 = xcx[12].r ; f4 = xcx[12].i ;
1422 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
1423 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
1424
1425 f1 = xcx[15].r ; f3 = xcx[15].i ; /* cos=1 sin=0 */
1426 f2 = xcx[14].r ; f4 = xcx[14].i ;
1427 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
1428 xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ;
1429
1430 f1 = xcx[17].r ; f3 = xcx[17].i ; /* cos=1 sin=0 */
1431 f2 = xcx[16].r ; f4 = xcx[16].i ;
1432 xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ;
1433 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
1434
1435 f1 = xcx[19].r ; f3 = xcx[19].i ; /* cos=1 sin=0 */
1436 f2 = xcx[18].r ; f4 = xcx[18].i ;
1437 xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ;
1438 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
1439
1440 f1 = xcx[21].r ; f3 = xcx[21].i ; /* cos=1 sin=0 */
1441 f2 = xcx[20].r ; f4 = xcx[20].i ;
1442 xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ;
1443 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
1444
1445 f1 = xcx[23].r ; f3 = xcx[23].i ; /* cos=1 sin=0 */
1446 f2 = xcx[22].r ; f4 = xcx[22].i ;
1447 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
1448 xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ;
1449
1450 f1 = xcx[25].r ; f3 = xcx[25].i ; /* cos=1 sin=0 */
1451 f2 = xcx[24].r ; f4 = xcx[24].i ;
1452 xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ;
1453 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
1454
1455 f1 = xcx[27].r ; f3 = xcx[27].i ; /* cos=1 sin=0 */
1456 f2 = xcx[26].r ; f4 = xcx[26].i ;
1457 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
1458 xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ;
1459
1460 f1 = xcx[29].r ; f3 = xcx[29].i ; /* cos=1 sin=0 */
1461 f2 = xcx[28].r ; f4 = xcx[28].i ;
1462 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
1463 xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ;
1464
1465 f1 = xcx[31].r ; f3 = xcx[31].i ; /* cos=1 sin=0 */
1466 f2 = xcx[30].r ; f4 = xcx[30].i ;
1467 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
1468 xcx[30].r = f2+f1 ; xcx[30].i = f4+f3 ;
1469
1470 f1 = xcx[33].r ; f3 = xcx[33].i ; /* cos=1 sin=0 */
1471 f2 = xcx[32].r ; f4 = xcx[32].i ;
1472 xcx[33].r = f2-f1 ; xcx[33].i = f4-f3 ;
1473 xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ;
1474
1475 f1 = xcx[35].r ; f3 = xcx[35].i ; /* cos=1 sin=0 */
1476 f2 = xcx[34].r ; f4 = xcx[34].i ;
1477 xcx[35].r = f2-f1 ; xcx[35].i = f4-f3 ;
1478 xcx[34].r = f2+f1 ; xcx[34].i = f4+f3 ;
1479
1480 f1 = xcx[37].r ; f3 = xcx[37].i ; /* cos=1 sin=0 */
1481 f2 = xcx[36].r ; f4 = xcx[36].i ;
1482 xcx[37].r = f2-f1 ; xcx[37].i = f4-f3 ;
1483 xcx[36].r = f2+f1 ; xcx[36].i = f4+f3 ;
1484
1485 f1 = xcx[39].r ; f3 = xcx[39].i ; /* cos=1 sin=0 */
1486 f2 = xcx[38].r ; f4 = xcx[38].i ;
1487 xcx[39].r = f2-f1 ; xcx[39].i = f4-f3 ;
1488 xcx[38].r = f2+f1 ; xcx[38].i = f4+f3 ;
1489
1490 f1 = xcx[41].r ; f3 = xcx[41].i ; /* cos=1 sin=0 */
1491 f2 = xcx[40].r ; f4 = xcx[40].i ;
1492 xcx[41].r = f2-f1 ; xcx[41].i = f4-f3 ;
1493 xcx[40].r = f2+f1 ; xcx[40].i = f4+f3 ;
1494
1495 f1 = xcx[43].r ; f3 = xcx[43].i ; /* cos=1 sin=0 */
1496 f2 = xcx[42].r ; f4 = xcx[42].i ;
1497 xcx[43].r = f2-f1 ; xcx[43].i = f4-f3 ;
1498 xcx[42].r = f2+f1 ; xcx[42].i = f4+f3 ;
1499
1500 f1 = xcx[45].r ; f3 = xcx[45].i ; /* cos=1 sin=0 */
1501 f2 = xcx[44].r ; f4 = xcx[44].i ;
1502 xcx[45].r = f2-f1 ; xcx[45].i = f4-f3 ;
1503 xcx[44].r = f2+f1 ; xcx[44].i = f4+f3 ;
1504
1505 f1 = xcx[47].r ; f3 = xcx[47].i ; /* cos=1 sin=0 */
1506 f2 = xcx[46].r ; f4 = xcx[46].i ;
1507 xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ;
1508 xcx[46].r = f2+f1 ; xcx[46].i = f4+f3 ;
1509
1510 f1 = xcx[49].r ; f3 = xcx[49].i ; /* cos=1 sin=0 */
1511 f2 = xcx[48].r ; f4 = xcx[48].i ;
1512 xcx[49].r = f2-f1 ; xcx[49].i = f4-f3 ;
1513 xcx[48].r = f2+f1 ; xcx[48].i = f4+f3 ;
1514
1515 f1 = xcx[51].r ; f3 = xcx[51].i ; /* cos=1 sin=0 */
1516 f2 = xcx[50].r ; f4 = xcx[50].i ;
1517 xcx[51].r = f2-f1 ; xcx[51].i = f4-f3 ;
1518 xcx[50].r = f2+f1 ; xcx[50].i = f4+f3 ;
1519
1520 f1 = xcx[53].r ; f3 = xcx[53].i ; /* cos=1 sin=0 */
1521 f2 = xcx[52].r ; f4 = xcx[52].i ;
1522 xcx[53].r = f2-f1 ; xcx[53].i = f4-f3 ;
1523 xcx[52].r = f2+f1 ; xcx[52].i = f4+f3 ;
1524
1525 f1 = xcx[55].r ; f3 = xcx[55].i ; /* cos=1 sin=0 */
1526 f2 = xcx[54].r ; f4 = xcx[54].i ;
1527 xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ;
1528 xcx[54].r = f2+f1 ; xcx[54].i = f4+f3 ;
1529
1530 f1 = xcx[57].r ; f3 = xcx[57].i ; /* cos=1 sin=0 */
1531 f2 = xcx[56].r ; f4 = xcx[56].i ;
1532 xcx[57].r = f2-f1 ; xcx[57].i = f4-f3 ;
1533 xcx[56].r = f2+f1 ; xcx[56].i = f4+f3 ;
1534
1535 f1 = xcx[59].r ; f3 = xcx[59].i ; /* cos=1 sin=0 */
1536 f2 = xcx[58].r ; f4 = xcx[58].i ;
1537 xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ;
1538 xcx[58].r = f2+f1 ; xcx[58].i = f4+f3 ;
1539
1540 f1 = xcx[61].r ; f3 = xcx[61].i ; /* cos=1 sin=0 */
1541 f2 = xcx[60].r ; f4 = xcx[60].i ;
1542 xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ;
1543 xcx[60].r = f2+f1 ; xcx[60].i = f4+f3 ;
1544
1545 f1 = xcx[63].r ; f3 = xcx[63].i ; /* cos=1 sin=0 */
1546 f2 = xcx[62].r ; f4 = xcx[62].i ;
1547 xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ;
1548 xcx[62].r = f2+f1 ; xcx[62].i = f4+f3 ;
1549
1550 f1 = xcx[2].r ; f3 = xcx[2].i ; /* cos=1 sin=0 */
1551 f2 = xcx[0].r ; f4 = xcx[0].i ;
1552 xcx[2].r = f2-f1 ; xcx[2].i = f4-f3 ;
1553 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
1554
1555 f1 = xcx[6].r ; f3 = xcx[6].i ; /* cos=1 sin=0 */
1556 f2 = xcx[4].r ; f4 = xcx[4].i ;
1557 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
1558 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
1559
1560 f1 = xcx[10].r ; f3 = xcx[10].i ; /* cos=1 sin=0 */
1561 f2 = xcx[8].r ; f4 = xcx[8].i ;
1562 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
1563 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
1564
1565 f1 = xcx[14].r ; f3 = xcx[14].i ; /* cos=1 sin=0 */
1566 f2 = xcx[12].r ; f4 = xcx[12].i ;
1567 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
1568 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
1569
1570 f1 = xcx[18].r ; f3 = xcx[18].i ; /* cos=1 sin=0 */
1571 f2 = xcx[16].r ; f4 = xcx[16].i ;
1572 xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ;
1573 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
1574
1575 f1 = xcx[22].r ; f3 = xcx[22].i ; /* cos=1 sin=0 */
1576 f2 = xcx[20].r ; f4 = xcx[20].i ;
1577 xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ;
1578 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
1579
1580 f1 = xcx[26].r ; f3 = xcx[26].i ; /* cos=1 sin=0 */
1581 f2 = xcx[24].r ; f4 = xcx[24].i ;
1582 xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ;
1583 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
1584
1585 f1 = xcx[30].r ; f3 = xcx[30].i ; /* cos=1 sin=0 */
1586 f2 = xcx[28].r ; f4 = xcx[28].i ;
1587 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
1588 xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ;
1589
1590 f1 = xcx[34].r ; f3 = xcx[34].i ; /* cos=1 sin=0 */
1591 f2 = xcx[32].r ; f4 = xcx[32].i ;
1592 xcx[34].r = f2-f1 ; xcx[34].i = f4-f3 ;
1593 xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ;
1594
1595 f1 = xcx[38].r ; f3 = xcx[38].i ; /* cos=1 sin=0 */
1596 f2 = xcx[36].r ; f4 = xcx[36].i ;
1597 xcx[38].r = f2-f1 ; xcx[38].i = f4-f3 ;
1598 xcx[36].r = f2+f1 ; xcx[36].i = f4+f3 ;
1599
1600 f1 = xcx[42].r ; f3 = xcx[42].i ; /* cos=1 sin=0 */
1601 f2 = xcx[40].r ; f4 = xcx[40].i ;
1602 xcx[42].r = f2-f1 ; xcx[42].i = f4-f3 ;
1603 xcx[40].r = f2+f1 ; xcx[40].i = f4+f3 ;
1604
1605 f1 = xcx[46].r ; f3 = xcx[46].i ; /* cos=1 sin=0 */
1606 f2 = xcx[44].r ; f4 = xcx[44].i ;
1607 xcx[46].r = f2-f1 ; xcx[46].i = f4-f3 ;
1608 xcx[44].r = f2+f1 ; xcx[44].i = f4+f3 ;
1609
1610 f1 = xcx[50].r ; f3 = xcx[50].i ; /* cos=1 sin=0 */
1611 f2 = xcx[48].r ; f4 = xcx[48].i ;
1612 xcx[50].r = f2-f1 ; xcx[50].i = f4-f3 ;
1613 xcx[48].r = f2+f1 ; xcx[48].i = f4+f3 ;
1614
1615 f1 = xcx[54].r ; f3 = xcx[54].i ; /* cos=1 sin=0 */
1616 f2 = xcx[52].r ; f4 = xcx[52].i ;
1617 xcx[54].r = f2-f1 ; xcx[54].i = f4-f3 ;
1618 xcx[52].r = f2+f1 ; xcx[52].i = f4+f3 ;
1619
1620 f1 = xcx[58].r ; f3 = xcx[58].i ; /* cos=1 sin=0 */
1621 f2 = xcx[56].r ; f4 = xcx[56].i ;
1622 xcx[58].r = f2-f1 ; xcx[58].i = f4-f3 ;
1623 xcx[56].r = f2+f1 ; xcx[56].i = f4+f3 ;
1624
1625 f1 = xcx[62].r ; f3 = xcx[62].i ; /* cos=1 sin=0 */
1626 f2 = xcx[60].r ; f4 = xcx[60].i ;
1627 xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ;
1628 xcx[60].r = f2+f1 ; xcx[60].i = f4+f3 ;
1629
1630 f1 = - xcx[3].i * csp[2].i ; /* cos=0 twiddles */
1631 f3 = xcx[3].r * csp[2].i ;
1632 f2 = xcx[1].r ; f4 = xcx[1].i ;
1633 xcx[3].r = f2-f1 ; xcx[3].i = f4-f3 ;
1634 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
1635
1636 f1 = - xcx[7].i * csp[2].i ; /* cos=0 twiddles */
1637 f3 = xcx[7].r * csp[2].i ;
1638 f2 = xcx[5].r ; f4 = xcx[5].i ;
1639 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
1640 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
1641
1642 f1 = - xcx[11].i * csp[2].i ; /* cos=0 twiddles */
1643 f3 = xcx[11].r * csp[2].i ;
1644 f2 = xcx[9].r ; f4 = xcx[9].i ;
1645 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
1646 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
1647
1648 f1 = - xcx[15].i * csp[2].i ; /* cos=0 twiddles */
1649 f3 = xcx[15].r * csp[2].i ;
1650 f2 = xcx[13].r ; f4 = xcx[13].i ;
1651 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
1652 xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ;
1653
1654 f1 = - xcx[19].i * csp[2].i ; /* cos=0 twiddles */
1655 f3 = xcx[19].r * csp[2].i ;
1656 f2 = xcx[17].r ; f4 = xcx[17].i ;
1657 xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ;
1658 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
1659
1660 f1 = - xcx[23].i * csp[2].i ; /* cos=0 twiddles */
1661 f3 = xcx[23].r * csp[2].i ;
1662 f2 = xcx[21].r ; f4 = xcx[21].i ;
1663 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
1664 xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ;
1665
1666 f1 = - xcx[27].i * csp[2].i ; /* cos=0 twiddles */
1667 f3 = xcx[27].r * csp[2].i ;
1668 f2 = xcx[25].r ; f4 = xcx[25].i ;
1669 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
1670 xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ;
1671
1672 f1 = - xcx[31].i * csp[2].i ; /* cos=0 twiddles */
1673 f3 = xcx[31].r * csp[2].i ;
1674 f2 = xcx[29].r ; f4 = xcx[29].i ;
1675 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
1676 xcx[29].r = f2+f1 ; xcx[29].i = f4+f3 ;
1677
1678 f1 = - xcx[35].i * csp[2].i ; /* cos=0 twiddles */
1679 f3 = xcx[35].r * csp[2].i ;
1680 f2 = xcx[33].r ; f4 = xcx[33].i ;
1681 xcx[35].r = f2-f1 ; xcx[35].i = f4-f3 ;
1682 xcx[33].r = f2+f1 ; xcx[33].i = f4+f3 ;
1683
1684 f1 = - xcx[39].i * csp[2].i ; /* cos=0 twiddles */
1685 f3 = xcx[39].r * csp[2].i ;
1686 f2 = xcx[37].r ; f4 = xcx[37].i ;
1687 xcx[39].r = f2-f1 ; xcx[39].i = f4-f3 ;
1688 xcx[37].r = f2+f1 ; xcx[37].i = f4+f3 ;
1689
1690 f1 = - xcx[43].i * csp[2].i ; /* cos=0 twiddles */
1691 f3 = xcx[43].r * csp[2].i ;
1692 f2 = xcx[41].r ; f4 = xcx[41].i ;
1693 xcx[43].r = f2-f1 ; xcx[43].i = f4-f3 ;
1694 xcx[41].r = f2+f1 ; xcx[41].i = f4+f3 ;
1695
1696 f1 = - xcx[47].i * csp[2].i ; /* cos=0 twiddles */
1697 f3 = xcx[47].r * csp[2].i ;
1698 f2 = xcx[45].r ; f4 = xcx[45].i ;
1699 xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ;
1700 xcx[45].r = f2+f1 ; xcx[45].i = f4+f3 ;
1701
1702 f1 = - xcx[51].i * csp[2].i ; /* cos=0 twiddles */
1703 f3 = xcx[51].r * csp[2].i ;
1704 f2 = xcx[49].r ; f4 = xcx[49].i ;
1705 xcx[51].r = f2-f1 ; xcx[51].i = f4-f3 ;
1706 xcx[49].r = f2+f1 ; xcx[49].i = f4+f3 ;
1707
1708 f1 = - xcx[55].i * csp[2].i ; /* cos=0 twiddles */
1709 f3 = xcx[55].r * csp[2].i ;
1710 f2 = xcx[53].r ; f4 = xcx[53].i ;
1711 xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ;
1712 xcx[53].r = f2+f1 ; xcx[53].i = f4+f3 ;
1713
1714 f1 = - xcx[59].i * csp[2].i ; /* cos=0 twiddles */
1715 f3 = xcx[59].r * csp[2].i ;
1716 f2 = xcx[57].r ; f4 = xcx[57].i ;
1717 xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ;
1718 xcx[57].r = f2+f1 ; xcx[57].i = f4+f3 ;
1719
1720 f1 = - xcx[63].i * csp[2].i ; /* cos=0 twiddles */
1721 f3 = xcx[63].r * csp[2].i ;
1722 f2 = xcx[61].r ; f4 = xcx[61].i ;
1723 xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ;
1724 xcx[61].r = f2+f1 ; xcx[61].i = f4+f3 ;
1725
1726 f1 = xcx[4].r ; f3 = xcx[4].i ; /* cos=1 sin=0 */
1727 f2 = xcx[0].r ; f4 = xcx[0].i ;
1728 xcx[4].r = f2-f1 ; xcx[4].i = f4-f3 ;
1729 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
1730
1731 f1 = xcx[12].r ; f3 = xcx[12].i ; /* cos=1 sin=0 */
1732 f2 = xcx[8].r ; f4 = xcx[8].i ;
1733 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
1734 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
1735
1736 f1 = xcx[20].r ; f3 = xcx[20].i ; /* cos=1 sin=0 */
1737 f2 = xcx[16].r ; f4 = xcx[16].i ;
1738 xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ;
1739 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
1740
1741 f1 = xcx[28].r ; f3 = xcx[28].i ; /* cos=1 sin=0 */
1742 f2 = xcx[24].r ; f4 = xcx[24].i ;
1743 xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ;
1744 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
1745
1746 f1 = xcx[36].r ; f3 = xcx[36].i ; /* cos=1 sin=0 */
1747 f2 = xcx[32].r ; f4 = xcx[32].i ;
1748 xcx[36].r = f2-f1 ; xcx[36].i = f4-f3 ;
1749 xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ;
1750
1751 f1 = xcx[44].r ; f3 = xcx[44].i ; /* cos=1 sin=0 */
1752 f2 = xcx[40].r ; f4 = xcx[40].i ;
1753 xcx[44].r = f2-f1 ; xcx[44].i = f4-f3 ;
1754 xcx[40].r = f2+f1 ; xcx[40].i = f4+f3 ;
1755
1756 f1 = xcx[52].r ; f3 = xcx[52].i ; /* cos=1 sin=0 */
1757 f2 = xcx[48].r ; f4 = xcx[48].i ;
1758 xcx[52].r = f2-f1 ; xcx[52].i = f4-f3 ;
1759 xcx[48].r = f2+f1 ; xcx[48].i = f4+f3 ;
1760
1761 f1 = xcx[60].r ; f3 = xcx[60].i ; /* cos=1 sin=0 */
1762 f2 = xcx[56].r ; f4 = xcx[56].i ;
1763 xcx[60].r = f2-f1 ; xcx[60].i = f4-f3 ;
1764 xcx[56].r = f2+f1 ; xcx[56].i = f4+f3 ;
1765
1766 f1 = xcx[5].r * csp[4].r - xcx[5].i * csp[4].i ; /* twiddles */
1767 f3 = xcx[5].r * csp[4].i + xcx[5].i * csp[4].r ;
1768 f2 = xcx[1].r ; f4 = xcx[1].i ;
1769 xcx[5].r = f2-f1 ; xcx[5].i = f4-f3 ;
1770 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
1771
1772 f1 = xcx[13].r * csp[4].r - xcx[13].i * csp[4].i ; /* twiddles */
1773 f3 = xcx[13].r * csp[4].i + xcx[13].i * csp[4].r ;
1774 f2 = xcx[9].r ; f4 = xcx[9].i ;
1775 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
1776 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
1777
1778 f1 = xcx[21].r * csp[4].r - xcx[21].i * csp[4].i ; /* twiddles */
1779 f3 = xcx[21].r * csp[4].i + xcx[21].i * csp[4].r ;
1780 f2 = xcx[17].r ; f4 = xcx[17].i ;
1781 xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ;
1782 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
1783
1784 f1 = xcx[29].r * csp[4].r - xcx[29].i * csp[4].i ; /* twiddles */
1785 f3 = xcx[29].r * csp[4].i + xcx[29].i * csp[4].r ;
1786 f2 = xcx[25].r ; f4 = xcx[25].i ;
1787 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
1788 xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ;
1789
1790 f1 = xcx[37].r * csp[4].r - xcx[37].i * csp[4].i ; /* twiddles */
1791 f3 = xcx[37].r * csp[4].i + xcx[37].i * csp[4].r ;
1792 f2 = xcx[33].r ; f4 = xcx[33].i ;
1793 xcx[37].r = f2-f1 ; xcx[37].i = f4-f3 ;
1794 xcx[33].r = f2+f1 ; xcx[33].i = f4+f3 ;
1795
1796 f1 = xcx[45].r * csp[4].r - xcx[45].i * csp[4].i ; /* twiddles */
1797 f3 = xcx[45].r * csp[4].i + xcx[45].i * csp[4].r ;
1798 f2 = xcx[41].r ; f4 = xcx[41].i ;
1799 xcx[45].r = f2-f1 ; xcx[45].i = f4-f3 ;
1800 xcx[41].r = f2+f1 ; xcx[41].i = f4+f3 ;
1801
1802 f1 = xcx[53].r * csp[4].r - xcx[53].i * csp[4].i ; /* twiddles */
1803 f3 = xcx[53].r * csp[4].i + xcx[53].i * csp[4].r ;
1804 f2 = xcx[49].r ; f4 = xcx[49].i ;
1805 xcx[53].r = f2-f1 ; xcx[53].i = f4-f3 ;
1806 xcx[49].r = f2+f1 ; xcx[49].i = f4+f3 ;
1807
1808 f1 = xcx[61].r * csp[4].r - xcx[61].i * csp[4].i ; /* twiddles */
1809 f3 = xcx[61].r * csp[4].i + xcx[61].i * csp[4].r ;
1810 f2 = xcx[57].r ; f4 = xcx[57].i ;
1811 xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ;
1812 xcx[57].r = f2+f1 ; xcx[57].i = f4+f3 ;
1813
1814 f1 = - xcx[6].i * csp[5].i ; /* cos=0 twiddles */
1815 f3 = xcx[6].r * csp[5].i ;
1816 f2 = xcx[2].r ; f4 = xcx[2].i ;
1817 xcx[6].r = f2-f1 ; xcx[6].i = f4-f3 ;
1818 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
1819
1820 f1 = - xcx[14].i * csp[5].i ; /* cos=0 twiddles */
1821 f3 = xcx[14].r * csp[5].i ;
1822 f2 = xcx[10].r ; f4 = xcx[10].i ;
1823 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
1824 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
1825
1826 f1 = - xcx[22].i * csp[5].i ; /* cos=0 twiddles */
1827 f3 = xcx[22].r * csp[5].i ;
1828 f2 = xcx[18].r ; f4 = xcx[18].i ;
1829 xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ;
1830 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
1831
1832 f1 = - xcx[30].i * csp[5].i ; /* cos=0 twiddles */
1833 f3 = xcx[30].r * csp[5].i ;
1834 f2 = xcx[26].r ; f4 = xcx[26].i ;
1835 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
1836 xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ;
1837
1838 f1 = - xcx[38].i * csp[5].i ; /* cos=0 twiddles */
1839 f3 = xcx[38].r * csp[5].i ;
1840 f2 = xcx[34].r ; f4 = xcx[34].i ;
1841 xcx[38].r = f2-f1 ; xcx[38].i = f4-f3 ;
1842 xcx[34].r = f2+f1 ; xcx[34].i = f4+f3 ;
1843
1844 f1 = - xcx[46].i * csp[5].i ; /* cos=0 twiddles */
1845 f3 = xcx[46].r * csp[5].i ;
1846 f2 = xcx[42].r ; f4 = xcx[42].i ;
1847 xcx[46].r = f2-f1 ; xcx[46].i = f4-f3 ;
1848 xcx[42].r = f2+f1 ; xcx[42].i = f4+f3 ;
1849
1850 f1 = - xcx[54].i * csp[5].i ; /* cos=0 twiddles */
1851 f3 = xcx[54].r * csp[5].i ;
1852 f2 = xcx[50].r ; f4 = xcx[50].i ;
1853 xcx[54].r = f2-f1 ; xcx[54].i = f4-f3 ;
1854 xcx[50].r = f2+f1 ; xcx[50].i = f4+f3 ;
1855
1856 f1 = - xcx[62].i * csp[5].i ; /* cos=0 twiddles */
1857 f3 = xcx[62].r * csp[5].i ;
1858 f2 = xcx[58].r ; f4 = xcx[58].i ;
1859 xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ;
1860 xcx[58].r = f2+f1 ; xcx[58].i = f4+f3 ;
1861
1862 f1 = xcx[7].r * csp[6].r - xcx[7].i * csp[6].i ; /* twiddles */
1863 f3 = xcx[7].r * csp[6].i + xcx[7].i * csp[6].r ;
1864 f2 = xcx[3].r ; f4 = xcx[3].i ;
1865 xcx[7].r = f2-f1 ; xcx[7].i = f4-f3 ;
1866 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
1867
1868 f1 = xcx[15].r * csp[6].r - xcx[15].i * csp[6].i ; /* twiddles */
1869 f3 = xcx[15].r * csp[6].i + xcx[15].i * csp[6].r ;
1870 f2 = xcx[11].r ; f4 = xcx[11].i ;
1871 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
1872 xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ;
1873
1874 f1 = xcx[23].r * csp[6].r - xcx[23].i * csp[6].i ; /* twiddles */
1875 f3 = xcx[23].r * csp[6].i + xcx[23].i * csp[6].r ;
1876 f2 = xcx[19].r ; f4 = xcx[19].i ;
1877 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
1878 xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ;
1879
1880 f1 = xcx[31].r * csp[6].r - xcx[31].i * csp[6].i ; /* twiddles */
1881 f3 = xcx[31].r * csp[6].i + xcx[31].i * csp[6].r ;
1882 f2 = xcx[27].r ; f4 = xcx[27].i ;
1883 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
1884 xcx[27].r = f2+f1 ; xcx[27].i = f4+f3 ;
1885
1886 f1 = xcx[39].r * csp[6].r - xcx[39].i * csp[6].i ; /* twiddles */
1887 f3 = xcx[39].r * csp[6].i + xcx[39].i * csp[6].r ;
1888 f2 = xcx[35].r ; f4 = xcx[35].i ;
1889 xcx[39].r = f2-f1 ; xcx[39].i = f4-f3 ;
1890 xcx[35].r = f2+f1 ; xcx[35].i = f4+f3 ;
1891
1892 f1 = xcx[47].r * csp[6].r - xcx[47].i * csp[6].i ; /* twiddles */
1893 f3 = xcx[47].r * csp[6].i + xcx[47].i * csp[6].r ;
1894 f2 = xcx[43].r ; f4 = xcx[43].i ;
1895 xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ;
1896 xcx[43].r = f2+f1 ; xcx[43].i = f4+f3 ;
1897
1898 f1 = xcx[55].r * csp[6].r - xcx[55].i * csp[6].i ; /* twiddles */
1899 f3 = xcx[55].r * csp[6].i + xcx[55].i * csp[6].r ;
1900 f2 = xcx[51].r ; f4 = xcx[51].i ;
1901 xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ;
1902 xcx[51].r = f2+f1 ; xcx[51].i = f4+f3 ;
1903
1904 f1 = xcx[63].r * csp[6].r - xcx[63].i * csp[6].i ; /* twiddles */
1905 f3 = xcx[63].r * csp[6].i + xcx[63].i * csp[6].r ;
1906 f2 = xcx[59].r ; f4 = xcx[59].i ;
1907 xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ;
1908 xcx[59].r = f2+f1 ; xcx[59].i = f4+f3 ;
1909
1910 f1 = xcx[8].r ; f3 = xcx[8].i ; /* cos=1 sin=0 */
1911 f2 = xcx[0].r ; f4 = xcx[0].i ;
1912 xcx[8].r = f2-f1 ; xcx[8].i = f4-f3 ;
1913 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
1914
1915 f1 = xcx[24].r ; f3 = xcx[24].i ; /* cos=1 sin=0 */
1916 f2 = xcx[16].r ; f4 = xcx[16].i ;
1917 xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ;
1918 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
1919
1920 f1 = xcx[40].r ; f3 = xcx[40].i ; /* cos=1 sin=0 */
1921 f2 = xcx[32].r ; f4 = xcx[32].i ;
1922 xcx[40].r = f2-f1 ; xcx[40].i = f4-f3 ;
1923 xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ;
1924
1925 f1 = xcx[56].r ; f3 = xcx[56].i ; /* cos=1 sin=0 */
1926 f2 = xcx[48].r ; f4 = xcx[48].i ;
1927 xcx[56].r = f2-f1 ; xcx[56].i = f4-f3 ;
1928 xcx[48].r = f2+f1 ; xcx[48].i = f4+f3 ;
1929
1930 f1 = xcx[9].r * csp[8].r - xcx[9].i * csp[8].i ; /* twiddles */
1931 f3 = xcx[9].r * csp[8].i + xcx[9].i * csp[8].r ;
1932 f2 = xcx[1].r ; f4 = xcx[1].i ;
1933 xcx[9].r = f2-f1 ; xcx[9].i = f4-f3 ;
1934 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
1935
1936 f1 = xcx[25].r * csp[8].r - xcx[25].i * csp[8].i ; /* twiddles */
1937 f3 = xcx[25].r * csp[8].i + xcx[25].i * csp[8].r ;
1938 f2 = xcx[17].r ; f4 = xcx[17].i ;
1939 xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ;
1940 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
1941
1942 f1 = xcx[41].r * csp[8].r - xcx[41].i * csp[8].i ; /* twiddles */
1943 f3 = xcx[41].r * csp[8].i + xcx[41].i * csp[8].r ;
1944 f2 = xcx[33].r ; f4 = xcx[33].i ;
1945 xcx[41].r = f2-f1 ; xcx[41].i = f4-f3 ;
1946 xcx[33].r = f2+f1 ; xcx[33].i = f4+f3 ;
1947
1948 f1 = xcx[57].r * csp[8].r - xcx[57].i * csp[8].i ; /* twiddles */
1949 f3 = xcx[57].r * csp[8].i + xcx[57].i * csp[8].r ;
1950 f2 = xcx[49].r ; f4 = xcx[49].i ;
1951 xcx[57].r = f2-f1 ; xcx[57].i = f4-f3 ;
1952 xcx[49].r = f2+f1 ; xcx[49].i = f4+f3 ;
1953
1954 f1 = xcx[10].r * csp[9].r - xcx[10].i * csp[9].i ; /* twiddles */
1955 f3 = xcx[10].r * csp[9].i + xcx[10].i * csp[9].r ;
1956 f2 = xcx[2].r ; f4 = xcx[2].i ;
1957 xcx[10].r = f2-f1 ; xcx[10].i = f4-f3 ;
1958 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
1959
1960 f1 = xcx[26].r * csp[9].r - xcx[26].i * csp[9].i ; /* twiddles */
1961 f3 = xcx[26].r * csp[9].i + xcx[26].i * csp[9].r ;
1962 f2 = xcx[18].r ; f4 = xcx[18].i ;
1963 xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ;
1964 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
1965
1966 f1 = xcx[42].r * csp[9].r - xcx[42].i * csp[9].i ; /* twiddles */
1967 f3 = xcx[42].r * csp[9].i + xcx[42].i * csp[9].r ;
1968 f2 = xcx[34].r ; f4 = xcx[34].i ;
1969 xcx[42].r = f2-f1 ; xcx[42].i = f4-f3 ;
1970 xcx[34].r = f2+f1 ; xcx[34].i = f4+f3 ;
1971
1972 f1 = xcx[58].r * csp[9].r - xcx[58].i * csp[9].i ; /* twiddles */
1973 f3 = xcx[58].r * csp[9].i + xcx[58].i * csp[9].r ;
1974 f2 = xcx[50].r ; f4 = xcx[50].i ;
1975 xcx[58].r = f2-f1 ; xcx[58].i = f4-f3 ;
1976 xcx[50].r = f2+f1 ; xcx[50].i = f4+f3 ;
1977
1978 f1 = xcx[11].r * csp[10].r - xcx[11].i * csp[10].i ; /* twiddles */
1979 f3 = xcx[11].r * csp[10].i + xcx[11].i * csp[10].r ;
1980 f2 = xcx[3].r ; f4 = xcx[3].i ;
1981 xcx[11].r = f2-f1 ; xcx[11].i = f4-f3 ;
1982 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
1983
1984 f1 = xcx[27].r * csp[10].r - xcx[27].i * csp[10].i ; /* twiddles */
1985 f3 = xcx[27].r * csp[10].i + xcx[27].i * csp[10].r ;
1986 f2 = xcx[19].r ; f4 = xcx[19].i ;
1987 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
1988 xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ;
1989
1990 f1 = xcx[43].r * csp[10].r - xcx[43].i * csp[10].i ; /* twiddles */
1991 f3 = xcx[43].r * csp[10].i + xcx[43].i * csp[10].r ;
1992 f2 = xcx[35].r ; f4 = xcx[35].i ;
1993 xcx[43].r = f2-f1 ; xcx[43].i = f4-f3 ;
1994 xcx[35].r = f2+f1 ; xcx[35].i = f4+f3 ;
1995
1996 f1 = xcx[59].r * csp[10].r - xcx[59].i * csp[10].i ; /* twiddles */
1997 f3 = xcx[59].r * csp[10].i + xcx[59].i * csp[10].r ;
1998 f2 = xcx[51].r ; f4 = xcx[51].i ;
1999 xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ;
2000 xcx[51].r = f2+f1 ; xcx[51].i = f4+f3 ;
2001
2002 f1 = - xcx[12].i * csp[11].i ; /* cos=0 twiddles */
2003 f3 = xcx[12].r * csp[11].i ;
2004 f2 = xcx[4].r ; f4 = xcx[4].i ;
2005 xcx[12].r = f2-f1 ; xcx[12].i = f4-f3 ;
2006 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
2007
2008 f1 = - xcx[28].i * csp[11].i ; /* cos=0 twiddles */
2009 f3 = xcx[28].r * csp[11].i ;
2010 f2 = xcx[20].r ; f4 = xcx[20].i ;
2011 xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ;
2012 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
2013
2014 f1 = - xcx[44].i * csp[11].i ; /* cos=0 twiddles */
2015 f3 = xcx[44].r * csp[11].i ;
2016 f2 = xcx[36].r ; f4 = xcx[36].i ;
2017 xcx[44].r = f2-f1 ; xcx[44].i = f4-f3 ;
2018 xcx[36].r = f2+f1 ; xcx[36].i = f4+f3 ;
2019
2020 f1 = - xcx[60].i * csp[11].i ; /* cos=0 twiddles */
2021 f3 = xcx[60].r * csp[11].i ;
2022 f2 = xcx[52].r ; f4 = xcx[52].i ;
2023 xcx[60].r = f2-f1 ; xcx[60].i = f4-f3 ;
2024 xcx[52].r = f2+f1 ; xcx[52].i = f4+f3 ;
2025
2026 f1 = xcx[13].r * csp[12].r - xcx[13].i * csp[12].i ; /* twiddles */
2027 f3 = xcx[13].r * csp[12].i + xcx[13].i * csp[12].r ;
2028 f2 = xcx[5].r ; f4 = xcx[5].i ;
2029 xcx[13].r = f2-f1 ; xcx[13].i = f4-f3 ;
2030 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
2031
2032 f1 = xcx[29].r * csp[12].r - xcx[29].i * csp[12].i ; /* twiddles */
2033 f3 = xcx[29].r * csp[12].i + xcx[29].i * csp[12].r ;
2034 f2 = xcx[21].r ; f4 = xcx[21].i ;
2035 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
2036 xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ;
2037
2038 f1 = xcx[45].r * csp[12].r - xcx[45].i * csp[12].i ; /* twiddles */
2039 f3 = xcx[45].r * csp[12].i + xcx[45].i * csp[12].r ;
2040 f2 = xcx[37].r ; f4 = xcx[37].i ;
2041 xcx[45].r = f2-f1 ; xcx[45].i = f4-f3 ;
2042 xcx[37].r = f2+f1 ; xcx[37].i = f4+f3 ;
2043
2044 f1 = xcx[61].r * csp[12].r - xcx[61].i * csp[12].i ; /* twiddles */
2045 f3 = xcx[61].r * csp[12].i + xcx[61].i * csp[12].r ;
2046 f2 = xcx[53].r ; f4 = xcx[53].i ;
2047 xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ;
2048 xcx[53].r = f2+f1 ; xcx[53].i = f4+f3 ;
2049
2050 f1 = xcx[14].r * csp[13].r - xcx[14].i * csp[13].i ; /* twiddles */
2051 f3 = xcx[14].r * csp[13].i + xcx[14].i * csp[13].r ;
2052 f2 = xcx[6].r ; f4 = xcx[6].i ;
2053 xcx[14].r = f2-f1 ; xcx[14].i = f4-f3 ;
2054 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
2055
2056 f1 = xcx[30].r * csp[13].r - xcx[30].i * csp[13].i ; /* twiddles */
2057 f3 = xcx[30].r * csp[13].i + xcx[30].i * csp[13].r ;
2058 f2 = xcx[22].r ; f4 = xcx[22].i ;
2059 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
2060 xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ;
2061
2062 f1 = xcx[46].r * csp[13].r - xcx[46].i * csp[13].i ; /* twiddles */
2063 f3 = xcx[46].r * csp[13].i + xcx[46].i * csp[13].r ;
2064 f2 = xcx[38].r ; f4 = xcx[38].i ;
2065 xcx[46].r = f2-f1 ; xcx[46].i = f4-f3 ;
2066 xcx[38].r = f2+f1 ; xcx[38].i = f4+f3 ;
2067
2068 f1 = xcx[62].r * csp[13].r - xcx[62].i * csp[13].i ; /* twiddles */
2069 f3 = xcx[62].r * csp[13].i + xcx[62].i * csp[13].r ;
2070 f2 = xcx[54].r ; f4 = xcx[54].i ;
2071 xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ;
2072 xcx[54].r = f2+f1 ; xcx[54].i = f4+f3 ;
2073
2074 f1 = xcx[15].r * csp[14].r - xcx[15].i * csp[14].i ; /* twiddles */
2075 f3 = xcx[15].r * csp[14].i + xcx[15].i * csp[14].r ;
2076 f2 = xcx[7].r ; f4 = xcx[7].i ;
2077 xcx[15].r = f2-f1 ; xcx[15].i = f4-f3 ;
2078 xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ;
2079
2080 f1 = xcx[31].r * csp[14].r - xcx[31].i * csp[14].i ; /* twiddles */
2081 f3 = xcx[31].r * csp[14].i + xcx[31].i * csp[14].r ;
2082 f2 = xcx[23].r ; f4 = xcx[23].i ;
2083 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
2084 xcx[23].r = f2+f1 ; xcx[23].i = f4+f3 ;
2085
2086 f1 = xcx[47].r * csp[14].r - xcx[47].i * csp[14].i ; /* twiddles */
2087 f3 = xcx[47].r * csp[14].i + xcx[47].i * csp[14].r ;
2088 f2 = xcx[39].r ; f4 = xcx[39].i ;
2089 xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ;
2090 xcx[39].r = f2+f1 ; xcx[39].i = f4+f3 ;
2091
2092 f1 = xcx[63].r * csp[14].r - xcx[63].i * csp[14].i ; /* twiddles */
2093 f3 = xcx[63].r * csp[14].i + xcx[63].i * csp[14].r ;
2094 f2 = xcx[55].r ; f4 = xcx[55].i ;
2095 xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ;
2096 xcx[55].r = f2+f1 ; xcx[55].i = f4+f3 ;
2097
2098 f1 = xcx[16].r ; f3 = xcx[16].i ; /* cos=1 sin=0 */
2099 f2 = xcx[0].r ; f4 = xcx[0].i ;
2100 xcx[16].r = f2-f1 ; xcx[16].i = f4-f3 ;
2101 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
2102
2103 f1 = xcx[48].r ; f3 = xcx[48].i ; /* cos=1 sin=0 */
2104 f2 = xcx[32].r ; f4 = xcx[32].i ;
2105 xcx[48].r = f2-f1 ; xcx[48].i = f4-f3 ;
2106 xcx[32].r = f2+f1 ; xcx[32].i = f4+f3 ;
2107
2108 f1 = xcx[17].r * csp[16].r - xcx[17].i * csp[16].i ; /* twiddles */
2109 f3 = xcx[17].r * csp[16].i + xcx[17].i * csp[16].r ;
2110 f2 = xcx[1].r ; f4 = xcx[1].i ;
2111 xcx[17].r = f2-f1 ; xcx[17].i = f4-f3 ;
2112 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
2113
2114 f1 = xcx[49].r * csp[16].r - xcx[49].i * csp[16].i ; /* twiddles */
2115 f3 = xcx[49].r * csp[16].i + xcx[49].i * csp[16].r ;
2116 f2 = xcx[33].r ; f4 = xcx[33].i ;
2117 xcx[49].r = f2-f1 ; xcx[49].i = f4-f3 ;
2118 xcx[33].r = f2+f1 ; xcx[33].i = f4+f3 ;
2119
2120 f1 = xcx[18].r * csp[17].r - xcx[18].i * csp[17].i ; /* twiddles */
2121 f3 = xcx[18].r * csp[17].i + xcx[18].i * csp[17].r ;
2122 f2 = xcx[2].r ; f4 = xcx[2].i ;
2123 xcx[18].r = f2-f1 ; xcx[18].i = f4-f3 ;
2124 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
2125
2126 f1 = xcx[50].r * csp[17].r - xcx[50].i * csp[17].i ; /* twiddles */
2127 f3 = xcx[50].r * csp[17].i + xcx[50].i * csp[17].r ;
2128 f2 = xcx[34].r ; f4 = xcx[34].i ;
2129 xcx[50].r = f2-f1 ; xcx[50].i = f4-f3 ;
2130 xcx[34].r = f2+f1 ; xcx[34].i = f4+f3 ;
2131
2132 f1 = xcx[19].r * csp[18].r - xcx[19].i * csp[18].i ; /* twiddles */
2133 f3 = xcx[19].r * csp[18].i + xcx[19].i * csp[18].r ;
2134 f2 = xcx[3].r ; f4 = xcx[3].i ;
2135 xcx[19].r = f2-f1 ; xcx[19].i = f4-f3 ;
2136 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
2137
2138 f1 = xcx[51].r * csp[18].r - xcx[51].i * csp[18].i ; /* twiddles */
2139 f3 = xcx[51].r * csp[18].i + xcx[51].i * csp[18].r ;
2140 f2 = xcx[35].r ; f4 = xcx[35].i ;
2141 xcx[51].r = f2-f1 ; xcx[51].i = f4-f3 ;
2142 xcx[35].r = f2+f1 ; xcx[35].i = f4+f3 ;
2143
2144 f1 = xcx[20].r * csp[19].r - xcx[20].i * csp[19].i ; /* twiddles */
2145 f3 = xcx[20].r * csp[19].i + xcx[20].i * csp[19].r ;
2146 f2 = xcx[4].r ; f4 = xcx[4].i ;
2147 xcx[20].r = f2-f1 ; xcx[20].i = f4-f3 ;
2148 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
2149
2150 f1 = xcx[52].r * csp[19].r - xcx[52].i * csp[19].i ; /* twiddles */
2151 f3 = xcx[52].r * csp[19].i + xcx[52].i * csp[19].r ;
2152 f2 = xcx[36].r ; f4 = xcx[36].i ;
2153 xcx[52].r = f2-f1 ; xcx[52].i = f4-f3 ;
2154 xcx[36].r = f2+f1 ; xcx[36].i = f4+f3 ;
2155
2156 f1 = xcx[21].r * csp[20].r - xcx[21].i * csp[20].i ; /* twiddles */
2157 f3 = xcx[21].r * csp[20].i + xcx[21].i * csp[20].r ;
2158 f2 = xcx[5].r ; f4 = xcx[5].i ;
2159 xcx[21].r = f2-f1 ; xcx[21].i = f4-f3 ;
2160 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
2161
2162 f1 = xcx[53].r * csp[20].r - xcx[53].i * csp[20].i ; /* twiddles */
2163 f3 = xcx[53].r * csp[20].i + xcx[53].i * csp[20].r ;
2164 f2 = xcx[37].r ; f4 = xcx[37].i ;
2165 xcx[53].r = f2-f1 ; xcx[53].i = f4-f3 ;
2166 xcx[37].r = f2+f1 ; xcx[37].i = f4+f3 ;
2167
2168 f1 = xcx[22].r * csp[21].r - xcx[22].i * csp[21].i ; /* twiddles */
2169 f3 = xcx[22].r * csp[21].i + xcx[22].i * csp[21].r ;
2170 f2 = xcx[6].r ; f4 = xcx[6].i ;
2171 xcx[22].r = f2-f1 ; xcx[22].i = f4-f3 ;
2172 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
2173
2174 f1 = xcx[54].r * csp[21].r - xcx[54].i * csp[21].i ; /* twiddles */
2175 f3 = xcx[54].r * csp[21].i + xcx[54].i * csp[21].r ;
2176 f2 = xcx[38].r ; f4 = xcx[38].i ;
2177 xcx[54].r = f2-f1 ; xcx[54].i = f4-f3 ;
2178 xcx[38].r = f2+f1 ; xcx[38].i = f4+f3 ;
2179
2180 f1 = xcx[23].r * csp[22].r - xcx[23].i * csp[22].i ; /* twiddles */
2181 f3 = xcx[23].r * csp[22].i + xcx[23].i * csp[22].r ;
2182 f2 = xcx[7].r ; f4 = xcx[7].i ;
2183 xcx[23].r = f2-f1 ; xcx[23].i = f4-f3 ;
2184 xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ;
2185
2186 f1 = xcx[55].r * csp[22].r - xcx[55].i * csp[22].i ; /* twiddles */
2187 f3 = xcx[55].r * csp[22].i + xcx[55].i * csp[22].r ;
2188 f2 = xcx[39].r ; f4 = xcx[39].i ;
2189 xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ;
2190 xcx[39].r = f2+f1 ; xcx[39].i = f4+f3 ;
2191
2192 f1 = - xcx[24].i * csp[23].i ; /* cos=0 twiddles */
2193 f3 = xcx[24].r * csp[23].i ;
2194 f2 = xcx[8].r ; f4 = xcx[8].i ;
2195 xcx[24].r = f2-f1 ; xcx[24].i = f4-f3 ;
2196 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
2197
2198 f1 = - xcx[56].i * csp[23].i ; /* cos=0 twiddles */
2199 f3 = xcx[56].r * csp[23].i ;
2200 f2 = xcx[40].r ; f4 = xcx[40].i ;
2201 xcx[56].r = f2-f1 ; xcx[56].i = f4-f3 ;
2202 xcx[40].r = f2+f1 ; xcx[40].i = f4+f3 ;
2203
2204 f1 = xcx[25].r * csp[24].r - xcx[25].i * csp[24].i ; /* twiddles */
2205 f3 = xcx[25].r * csp[24].i + xcx[25].i * csp[24].r ;
2206 f2 = xcx[9].r ; f4 = xcx[9].i ;
2207 xcx[25].r = f2-f1 ; xcx[25].i = f4-f3 ;
2208 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
2209
2210 f1 = xcx[57].r * csp[24].r - xcx[57].i * csp[24].i ; /* twiddles */
2211 f3 = xcx[57].r * csp[24].i + xcx[57].i * csp[24].r ;
2212 f2 = xcx[41].r ; f4 = xcx[41].i ;
2213 xcx[57].r = f2-f1 ; xcx[57].i = f4-f3 ;
2214 xcx[41].r = f2+f1 ; xcx[41].i = f4+f3 ;
2215
2216 f1 = xcx[26].r * csp[25].r - xcx[26].i * csp[25].i ; /* twiddles */
2217 f3 = xcx[26].r * csp[25].i + xcx[26].i * csp[25].r ;
2218 f2 = xcx[10].r ; f4 = xcx[10].i ;
2219 xcx[26].r = f2-f1 ; xcx[26].i = f4-f3 ;
2220 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
2221
2222 f1 = xcx[58].r * csp[25].r - xcx[58].i * csp[25].i ; /* twiddles */
2223 f3 = xcx[58].r * csp[25].i + xcx[58].i * csp[25].r ;
2224 f2 = xcx[42].r ; f4 = xcx[42].i ;
2225 xcx[58].r = f2-f1 ; xcx[58].i = f4-f3 ;
2226 xcx[42].r = f2+f1 ; xcx[42].i = f4+f3 ;
2227
2228 f1 = xcx[27].r * csp[26].r - xcx[27].i * csp[26].i ; /* twiddles */
2229 f3 = xcx[27].r * csp[26].i + xcx[27].i * csp[26].r ;
2230 f2 = xcx[11].r ; f4 = xcx[11].i ;
2231 xcx[27].r = f2-f1 ; xcx[27].i = f4-f3 ;
2232 xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ;
2233
2234 f1 = xcx[59].r * csp[26].r - xcx[59].i * csp[26].i ; /* twiddles */
2235 f3 = xcx[59].r * csp[26].i + xcx[59].i * csp[26].r ;
2236 f2 = xcx[43].r ; f4 = xcx[43].i ;
2237 xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ;
2238 xcx[43].r = f2+f1 ; xcx[43].i = f4+f3 ;
2239
2240 f1 = xcx[28].r * csp[27].r - xcx[28].i * csp[27].i ; /* twiddles */
2241 f3 = xcx[28].r * csp[27].i + xcx[28].i * csp[27].r ;
2242 f2 = xcx[12].r ; f4 = xcx[12].i ;
2243 xcx[28].r = f2-f1 ; xcx[28].i = f4-f3 ;
2244 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
2245
2246 f1 = xcx[60].r * csp[27].r - xcx[60].i * csp[27].i ; /* twiddles */
2247 f3 = xcx[60].r * csp[27].i + xcx[60].i * csp[27].r ;
2248 f2 = xcx[44].r ; f4 = xcx[44].i ;
2249 xcx[60].r = f2-f1 ; xcx[60].i = f4-f3 ;
2250 xcx[44].r = f2+f1 ; xcx[44].i = f4+f3 ;
2251
2252 f1 = xcx[29].r * csp[28].r - xcx[29].i * csp[28].i ; /* twiddles */
2253 f3 = xcx[29].r * csp[28].i + xcx[29].i * csp[28].r ;
2254 f2 = xcx[13].r ; f4 = xcx[13].i ;
2255 xcx[29].r = f2-f1 ; xcx[29].i = f4-f3 ;
2256 xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ;
2257
2258 f1 = xcx[61].r * csp[28].r - xcx[61].i * csp[28].i ; /* twiddles */
2259 f3 = xcx[61].r * csp[28].i + xcx[61].i * csp[28].r ;
2260 f2 = xcx[45].r ; f4 = xcx[45].i ;
2261 xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ;
2262 xcx[45].r = f2+f1 ; xcx[45].i = f4+f3 ;
2263
2264 f1 = xcx[30].r * csp[29].r - xcx[30].i * csp[29].i ; /* twiddles */
2265 f3 = xcx[30].r * csp[29].i + xcx[30].i * csp[29].r ;
2266 f2 = xcx[14].r ; f4 = xcx[14].i ;
2267 xcx[30].r = f2-f1 ; xcx[30].i = f4-f3 ;
2268 xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ;
2269
2270 f1 = xcx[62].r * csp[29].r - xcx[62].i * csp[29].i ; /* twiddles */
2271 f3 = xcx[62].r * csp[29].i + xcx[62].i * csp[29].r ;
2272 f2 = xcx[46].r ; f4 = xcx[46].i ;
2273 xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ;
2274 xcx[46].r = f2+f1 ; xcx[46].i = f4+f3 ;
2275
2276 f1 = xcx[31].r * csp[30].r - xcx[31].i * csp[30].i ; /* twiddles */
2277 f3 = xcx[31].r * csp[30].i + xcx[31].i * csp[30].r ;
2278 f2 = xcx[15].r ; f4 = xcx[15].i ;
2279 xcx[31].r = f2-f1 ; xcx[31].i = f4-f3 ;
2280 xcx[15].r = f2+f1 ; xcx[15].i = f4+f3 ;
2281
2282 f1 = xcx[63].r * csp[30].r - xcx[63].i * csp[30].i ; /* twiddles */
2283 f3 = xcx[63].r * csp[30].i + xcx[63].i * csp[30].r ;
2284 f2 = xcx[47].r ; f4 = xcx[47].i ;
2285 xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ;
2286 xcx[47].r = f2+f1 ; xcx[47].i = f4+f3 ;
2287
2288 f1 = xcx[32].r ; f3 = xcx[32].i ; /* cos=1 sin=0 */
2289 f2 = xcx[0].r ; f4 = xcx[0].i ;
2290 xcx[32].r = f2-f1 ; xcx[32].i = f4-f3 ;
2291 xcx[0].r = f2+f1 ; xcx[0].i = f4+f3 ;
2292
2293 f1 = xcx[33].r * csp[32].r - xcx[33].i * csp[32].i ; /* twiddles */
2294 f3 = xcx[33].r * csp[32].i + xcx[33].i * csp[32].r ;
2295 f2 = xcx[1].r ; f4 = xcx[1].i ;
2296 xcx[33].r = f2-f1 ; xcx[33].i = f4-f3 ;
2297 xcx[1].r = f2+f1 ; xcx[1].i = f4+f3 ;
2298
2299 f1 = xcx[34].r * csp[33].r - xcx[34].i * csp[33].i ; /* twiddles */
2300 f3 = xcx[34].r * csp[33].i + xcx[34].i * csp[33].r ;
2301 f2 = xcx[2].r ; f4 = xcx[2].i ;
2302 xcx[34].r = f2-f1 ; xcx[34].i = f4-f3 ;
2303 xcx[2].r = f2+f1 ; xcx[2].i = f4+f3 ;
2304
2305 f1 = xcx[35].r * csp[34].r - xcx[35].i * csp[34].i ; /* twiddles */
2306 f3 = xcx[35].r * csp[34].i + xcx[35].i * csp[34].r ;
2307 f2 = xcx[3].r ; f4 = xcx[3].i ;
2308 xcx[35].r = f2-f1 ; xcx[35].i = f4-f3 ;
2309 xcx[3].r = f2+f1 ; xcx[3].i = f4+f3 ;
2310
2311 f1 = xcx[36].r * csp[35].r - xcx[36].i * csp[35].i ; /* twiddles */
2312 f3 = xcx[36].r * csp[35].i + xcx[36].i * csp[35].r ;
2313 f2 = xcx[4].r ; f4 = xcx[4].i ;
2314 xcx[36].r = f2-f1 ; xcx[36].i = f4-f3 ;
2315 xcx[4].r = f2+f1 ; xcx[4].i = f4+f3 ;
2316
2317 f1 = xcx[37].r * csp[36].r - xcx[37].i * csp[36].i ; /* twiddles */
2318 f3 = xcx[37].r * csp[36].i + xcx[37].i * csp[36].r ;
2319 f2 = xcx[5].r ; f4 = xcx[5].i ;
2320 xcx[37].r = f2-f1 ; xcx[37].i = f4-f3 ;
2321 xcx[5].r = f2+f1 ; xcx[5].i = f4+f3 ;
2322
2323 f1 = xcx[38].r * csp[37].r - xcx[38].i * csp[37].i ; /* twiddles */
2324 f3 = xcx[38].r * csp[37].i + xcx[38].i * csp[37].r ;
2325 f2 = xcx[6].r ; f4 = xcx[6].i ;
2326 xcx[38].r = f2-f1 ; xcx[38].i = f4-f3 ;
2327 xcx[6].r = f2+f1 ; xcx[6].i = f4+f3 ;
2328
2329 f1 = xcx[39].r * csp[38].r - xcx[39].i * csp[38].i ; /* twiddles */
2330 f3 = xcx[39].r * csp[38].i + xcx[39].i * csp[38].r ;
2331 f2 = xcx[7].r ; f4 = xcx[7].i ;
2332 xcx[39].r = f2-f1 ; xcx[39].i = f4-f3 ;
2333 xcx[7].r = f2+f1 ; xcx[7].i = f4+f3 ;
2334
2335 f1 = xcx[40].r * csp[39].r - xcx[40].i * csp[39].i ; /* twiddles */
2336 f3 = xcx[40].r * csp[39].i + xcx[40].i * csp[39].r ;
2337 f2 = xcx[8].r ; f4 = xcx[8].i ;
2338 xcx[40].r = f2-f1 ; xcx[40].i = f4-f3 ;
2339 xcx[8].r = f2+f1 ; xcx[8].i = f4+f3 ;
2340
2341 f1 = xcx[41].r * csp[40].r - xcx[41].i * csp[40].i ; /* twiddles */
2342 f3 = xcx[41].r * csp[40].i + xcx[41].i * csp[40].r ;
2343 f2 = xcx[9].r ; f4 = xcx[9].i ;
2344 xcx[41].r = f2-f1 ; xcx[41].i = f4-f3 ;
2345 xcx[9].r = f2+f1 ; xcx[9].i = f4+f3 ;
2346
2347 f1 = xcx[42].r * csp[41].r - xcx[42].i * csp[41].i ; /* twiddles */
2348 f3 = xcx[42].r * csp[41].i + xcx[42].i * csp[41].r ;
2349 f2 = xcx[10].r ; f4 = xcx[10].i ;
2350 xcx[42].r = f2-f1 ; xcx[42].i = f4-f3 ;
2351 xcx[10].r = f2+f1 ; xcx[10].i = f4+f3 ;
2352
2353 f1 = xcx[43].r * csp[42].r - xcx[43].i * csp[42].i ; /* twiddles */
2354 f3 = xcx[43].r * csp[42].i + xcx[43].i * csp[42].r ;
2355 f2 = xcx[11].r ; f4 = xcx[11].i ;
2356 xcx[43].r = f2-f1 ; xcx[43].i = f4-f3 ;
2357 xcx[11].r = f2+f1 ; xcx[11].i = f4+f3 ;
2358
2359 f1 = xcx[44].r * csp[43].r - xcx[44].i * csp[43].i ; /* twiddles */
2360 f3 = xcx[44].r * csp[43].i + xcx[44].i * csp[43].r ;
2361 f2 = xcx[12].r ; f4 = xcx[12].i ;
2362 xcx[44].r = f2-f1 ; xcx[44].i = f4-f3 ;
2363 xcx[12].r = f2+f1 ; xcx[12].i = f4+f3 ;
2364
2365 f1 = xcx[45].r * csp[44].r - xcx[45].i * csp[44].i ; /* twiddles */
2366 f3 = xcx[45].r * csp[44].i + xcx[45].i * csp[44].r ;
2367 f2 = xcx[13].r ; f4 = xcx[13].i ;
2368 xcx[45].r = f2-f1 ; xcx[45].i = f4-f3 ;
2369 xcx[13].r = f2+f1 ; xcx[13].i = f4+f3 ;
2370
2371 f1 = xcx[46].r * csp[45].r - xcx[46].i * csp[45].i ; /* twiddles */
2372 f3 = xcx[46].r * csp[45].i + xcx[46].i * csp[45].r ;
2373 f2 = xcx[14].r ; f4 = xcx[14].i ;
2374 xcx[46].r = f2-f1 ; xcx[46].i = f4-f3 ;
2375 xcx[14].r = f2+f1 ; xcx[14].i = f4+f3 ;
2376
2377 f1 = xcx[47].r * csp[46].r - xcx[47].i * csp[46].i ; /* twiddles */
2378 f3 = xcx[47].r * csp[46].i + xcx[47].i * csp[46].r ;
2379 f2 = xcx[15].r ; f4 = xcx[15].i ;
2380 xcx[47].r = f2-f1 ; xcx[47].i = f4-f3 ;
2381 xcx[15].r = f2+f1 ; xcx[15].i = f4+f3 ;
2382
2383 f1 = - xcx[48].i * csp[47].i ; /* cos=0 twiddles */
2384 f3 = xcx[48].r * csp[47].i ;
2385 f2 = xcx[16].r ; f4 = xcx[16].i ;
2386 xcx[48].r = f2-f1 ; xcx[48].i = f4-f3 ;
2387 xcx[16].r = f2+f1 ; xcx[16].i = f4+f3 ;
2388
2389 f1 = xcx[49].r * csp[48].r - xcx[49].i * csp[48].i ; /* twiddles */
2390 f3 = xcx[49].r * csp[48].i + xcx[49].i * csp[48].r ;
2391 f2 = xcx[17].r ; f4 = xcx[17].i ;
2392 xcx[49].r = f2-f1 ; xcx[49].i = f4-f3 ;
2393 xcx[17].r = f2+f1 ; xcx[17].i = f4+f3 ;
2394
2395 f1 = xcx[50].r * csp[49].r - xcx[50].i * csp[49].i ; /* twiddles */
2396 f3 = xcx[50].r * csp[49].i + xcx[50].i * csp[49].r ;
2397 f2 = xcx[18].r ; f4 = xcx[18].i ;
2398 xcx[50].r = f2-f1 ; xcx[50].i = f4-f3 ;
2399 xcx[18].r = f2+f1 ; xcx[18].i = f4+f3 ;
2400
2401 f1 = xcx[51].r * csp[50].r - xcx[51].i * csp[50].i ; /* twiddles */
2402 f3 = xcx[51].r * csp[50].i + xcx[51].i * csp[50].r ;
2403 f2 = xcx[19].r ; f4 = xcx[19].i ;
2404 xcx[51].r = f2-f1 ; xcx[51].i = f4-f3 ;
2405 xcx[19].r = f2+f1 ; xcx[19].i = f4+f3 ;
2406
2407 f1 = xcx[52].r * csp[51].r - xcx[52].i * csp[51].i ; /* twiddles */
2408 f3 = xcx[52].r * csp[51].i + xcx[52].i * csp[51].r ;
2409 f2 = xcx[20].r ; f4 = xcx[20].i ;
2410 xcx[52].r = f2-f1 ; xcx[52].i = f4-f3 ;
2411 xcx[20].r = f2+f1 ; xcx[20].i = f4+f3 ;
2412
2413 f1 = xcx[53].r * csp[52].r - xcx[53].i * csp[52].i ; /* twiddles */
2414 f3 = xcx[53].r * csp[52].i + xcx[53].i * csp[52].r ;
2415 f2 = xcx[21].r ; f4 = xcx[21].i ;
2416 xcx[53].r = f2-f1 ; xcx[53].i = f4-f3 ;
2417 xcx[21].r = f2+f1 ; xcx[21].i = f4+f3 ;
2418
2419 f1 = xcx[54].r * csp[53].r - xcx[54].i * csp[53].i ; /* twiddles */
2420 f3 = xcx[54].r * csp[53].i + xcx[54].i * csp[53].r ;
2421 f2 = xcx[22].r ; f4 = xcx[22].i ;
2422 xcx[54].r = f2-f1 ; xcx[54].i = f4-f3 ;
2423 xcx[22].r = f2+f1 ; xcx[22].i = f4+f3 ;
2424
2425 f1 = xcx[55].r * csp[54].r - xcx[55].i * csp[54].i ; /* twiddles */
2426 f3 = xcx[55].r * csp[54].i + xcx[55].i * csp[54].r ;
2427 f2 = xcx[23].r ; f4 = xcx[23].i ;
2428 xcx[55].r = f2-f1 ; xcx[55].i = f4-f3 ;
2429 xcx[23].r = f2+f1 ; xcx[23].i = f4+f3 ;
2430
2431 f1 = xcx[56].r * csp[55].r - xcx[56].i * csp[55].i ; /* twiddles */
2432 f3 = xcx[56].r * csp[55].i + xcx[56].i * csp[55].r ;
2433 f2 = xcx[24].r ; f4 = xcx[24].i ;
2434 xcx[56].r = f2-f1 ; xcx[56].i = f4-f3 ;
2435 xcx[24].r = f2+f1 ; xcx[24].i = f4+f3 ;
2436
2437 f1 = xcx[57].r * csp[56].r - xcx[57].i * csp[56].i ; /* twiddles */
2438 f3 = xcx[57].r * csp[56].i + xcx[57].i * csp[56].r ;
2439 f2 = xcx[25].r ; f4 = xcx[25].i ;
2440 xcx[57].r = f2-f1 ; xcx[57].i = f4-f3 ;
2441 xcx[25].r = f2+f1 ; xcx[25].i = f4+f3 ;
2442
2443 f1 = xcx[58].r * csp[57].r - xcx[58].i * csp[57].i ; /* twiddles */
2444 f3 = xcx[58].r * csp[57].i + xcx[58].i * csp[57].r ;
2445 f2 = xcx[26].r ; f4 = xcx[26].i ;
2446 xcx[58].r = f2-f1 ; xcx[58].i = f4-f3 ;
2447 xcx[26].r = f2+f1 ; xcx[26].i = f4+f3 ;
2448
2449 f1 = xcx[59].r * csp[58].r - xcx[59].i * csp[58].i ; /* twiddles */
2450 f3 = xcx[59].r * csp[58].i + xcx[59].i * csp[58].r ;
2451 f2 = xcx[27].r ; f4 = xcx[27].i ;
2452 xcx[59].r = f2-f1 ; xcx[59].i = f4-f3 ;
2453 xcx[27].r = f2+f1 ; xcx[27].i = f4+f3 ;
2454
2455 f1 = xcx[60].r * csp[59].r - xcx[60].i * csp[59].i ; /* twiddles */
2456 f3 = xcx[60].r * csp[59].i + xcx[60].i * csp[59].r ;
2457 f2 = xcx[28].r ; f4 = xcx[28].i ;
2458 xcx[60].r = f2-f1 ; xcx[60].i = f4-f3 ;
2459 xcx[28].r = f2+f1 ; xcx[28].i = f4+f3 ;
2460
2461 f1 = xcx[61].r * csp[60].r - xcx[61].i * csp[60].i ; /* twiddles */
2462 f3 = xcx[61].r * csp[60].i + xcx[61].i * csp[60].r ;
2463 f2 = xcx[29].r ; f4 = xcx[29].i ;
2464 xcx[61].r = f2-f1 ; xcx[61].i = f4-f3 ;
2465 xcx[29].r = f2+f1 ; xcx[29].i = f4+f3 ;
2466
2467 f1 = xcx[62].r * csp[61].r - xcx[62].i * csp[61].i ; /* twiddles */
2468 f3 = xcx[62].r * csp[61].i + xcx[62].i * csp[61].r ;
2469 f2 = xcx[30].r ; f4 = xcx[30].i ;
2470 xcx[62].r = f2-f1 ; xcx[62].i = f4-f3 ;
2471 xcx[30].r = f2+f1 ; xcx[30].i = f4+f3 ;
2472
2473 f1 = xcx[63].r * csp[62].r - xcx[63].i * csp[62].i ; /* twiddles */
2474 f3 = xcx[63].r * csp[62].i + xcx[63].i * csp[62].r ;
2475 f2 = xcx[31].r ; f4 = xcx[31].i ;
2476 xcx[63].r = f2-f1 ; xcx[63].i = f4-f3 ;
2477 xcx[31].r = f2+f1 ; xcx[31].i = f4+f3 ;
2478
2479 return ;
2480 }
2481
2482 /*================================================================*/
2483
2484 /*------------------------------------------------------------------
2485 fft_4dec: do a decimation by 4.
2486 At most RMAX levels of recursion are allowed.
2487 --------------------------------------------------------------------*/
2488
fft_4dec_OMP(int mode,int idim,complex * xc)2489 static void fft_4dec_OMP( int mode , int idim , complex *xc )
2490 {
2491 DECLARE_ithr ;
2492 int rec , mold ;
2493 int N=idim , M=idim/4 , M2=2*M , M3=3*M ;
2494 complex *cs , *aa , *bb , *cc , *dd ;
2495 int k , i ;
2496 float aar,aai, tr,ti, bbr,bbi, ccr,cci , ddr,ddi , t1,t2 ,
2497 acpr,acmr , bdpr,bdmr , acpi,acmi , bdpi,bdmi ;
2498
2499 /* ININFO_message("fft_4dec_OMP(%d)",idim) ; */
2500
2501 /*-- initialize cosine/sine table and memory space --*/
2502
2503 rec = CFO_recX[ithr] ; /* recursion level (0 .. RMAX-1) */
2504 if( rec >= RMAX ){
2505 fprintf(stderr,"\n*** fft_4dec_OMP: too many recursions with idim=%d\n",idim); EXIT(1);
2506 }
2507
2508 CFO_assign(cs,M3) ; /* assign workspace of length M3 to cs pointer */
2509 CFO_assign(aa,M) ; /* at this recursion level for this thread (etc) */
2510 CFO_assign(bb,M) ;
2511 CFO_assign(cc,M) ;
2512 CFO_assign(dd,M) ;
2513
2514 mold = CFO_rmold4[rec][ithr] ; /* what was done before at this recursion level */
2515
2516 { if( M != mold ){
2517 double th = (2.0*PI/N) ;
2518 cs[0].r = 1.0 ; cs[0].i = 0.0 ;
2519 for( k=1 ; k < M3 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); }
2520 CFO_rmold4[rec][ithr] = M ;
2521 }}
2522
2523 /*-- load subarrays, and FFT each one --*/
2524
2525 for( i=k=0; k < M; k++ ){
2526 aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++];
2527 }
2528
2529 CFO_recX[ithr]++ ; /* bump recursion level */
2530 switch( N ){
2531
2532 default: fprintf(stderr,"\n*** Illegal call to fft_4dec: N=%d\n",N);
2533 EXIT(1) ;
2534
2535 case 128: fft32(mode,aa); fft32(mode,bb); /* no recursion */
2536 fft32(mode,cc); fft32(mode,dd); break;
2537
2538 case 256: fft64(mode,aa); fft64(mode,bb); /* no recursion */
2539 fft64(mode,cc); fft64(mode,dd); break;
2540
2541 case 512: fft128(mode,aa); fft128(mode,bb); /* recurse once */
2542 fft128(mode,cc); fft128(mode,dd); break;
2543
2544 case 1024: fft256(mode,aa); fft256(mode,bb); /* recurse once */
2545 fft256(mode,cc); fft256(mode,dd); break;
2546
2547 #if 1
2548 case 2048: fft512(mode,aa); fft512(mode,bb); /* recurse twice */
2549 fft512(mode,cc); fft512(mode,dd); break;
2550
2551 case 4096: fft1024(mode,aa); fft1024(mode,bb); /* recurse twice */
2552 fft1024(mode,cc); fft1024(mode,dd); break;
2553
2554 case 8192: fft2048(mode,aa); fft2048(mode,bb); /* recurse 3 times */
2555 fft2048(mode,cc); fft2048(mode,dd); break;
2556
2557 case 16384: fft4096(mode,aa); fft4096(mode,bb); /* recurse 3 times */
2558 fft4096(mode,cc); fft4096(mode,dd); break;
2559
2560 case 32768: fft8192(mode,aa); fft8192(mode,bb); /* recurse 4 times */
2561 fft8192(mode,cc); fft8192(mode,dd); break;
2562
2563 case 65536: fft16384(mode,aa); fft16384(mode,bb); /* recurse 4 times */
2564 fft16384(mode,cc); fft16384(mode,dd); break;
2565 #endif
2566 }
2567 CFO_recX[ithr]-- ; /* reduce recursion level */
2568
2569 /*-- recombination --*/
2570
2571 if( mode > 0 ){
2572 for( k=0 ; k < M ; k++ ){
2573 t1 = bb[k].r; t2 = bb[k].i;
2574 tr = cs[k].r; ti = cs[k].i;
2575 bbr = tr*t1 - ti*t2 ; bbi = tr*t2 + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */
2576
2577 t1 = cc[ k].r; t2 = cc[ k].i;
2578 tr = cs[2*k].r; ti = cs[2*k].i;
2579 ccr = tr*t1 - ti*t2 ; cci = tr*t2 + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */
2580
2581 t1 = dd[ k].r; t2 = dd[ k].i;
2582 tr = cs[3*k].r; ti = cs[3*k].i;
2583 ddr = tr*t1 - ti*t2 ; ddi = tr*t2 + ti*t1 ; /* d[k]*exp(+6*Pi*k/N) */
2584
2585 aar = aa[k].r; aai = aa[k].i; /* a[k] */
2586
2587 acpr = aar + ccr ; acmr = aar - ccr ;
2588 bdpr = bbr + ddr ; bdmr = bbr - ddr ;
2589 acpi = aai + cci ; acmi = aai - cci ;
2590 bdpi = bbi + ddi ; bdmi = bbi - ddi ;
2591
2592 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ;
2593 xc[k+M ].r = acmr-bdmi ; xc[k+M ].i = acmi+bdmr ;
2594 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
2595 xc[k+M3].r = acmr+bdmi ; xc[k+M3].i = acmi-bdmr ;
2596 }
2597 } else {
2598 for( k=0 ; k < M ; k++ ){
2599 t1 = bb[k].r; t2 = bb[k].i;
2600 tr = cs[k].r; ti = cs[k].i;
2601 bbr = tr*t1 + ti*t2 ; bbi = tr*t2 - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */
2602
2603 t1 = cc[ k].r; t2 = cc[ k].i;
2604 tr = cs[2*k].r; ti = cs[2*k].i;
2605 ccr = tr*t1 + ti*t2 ; cci = tr*t2 - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */
2606
2607 t1 = dd[ k].r; t2 = dd[ k].i;
2608 tr = cs[3*k].r; ti = cs[3*k].i;
2609 ddr = tr*t1 + ti*t2 ; ddi = tr*t2 - ti*t1 ; /* d[k]*exp(-6*Pi*k/N) */
2610
2611 aar = aa[k].r; aai = aa[k].i; /* a[k] */
2612
2613 acpr = aar + ccr ; acmr = aar - ccr ;
2614 bdpr = bbr + ddr ; bdmr = bbr - ddr ;
2615 acpi = aai + cci ; acmi = aai - cci ;
2616 bdpi = bbi + ddi ; bdmi = bbi - ddi ;
2617
2618 xc[k ].r = acpr+bdpr ; xc[k ].i = acpi+bdpi ;
2619 xc[k+M ].r = acmr+bdmi ; xc[k+M ].i = acmi-bdmr ;
2620 xc[k+M2].r = acpr-bdpr ; xc[k+M2].i = acpi-bdpi ;
2621 xc[k+M3].r = acmr-bdmi ; xc[k+M3].i = acmi+bdmr ;
2622 }
2623 }
2624
2625 return ;
2626 }
2627
2628 /*=======================================================================
2629 The following radix-3 and radix-5 routines are by RWCox -- Aug 1999.
2630 They are not inside unrolled code.
2631 =========================================================================*/
2632
2633 /*------------------------------------------------------------------
2634 fft_3dec: do a decimation-by-3, plus a power-of-2.
2635 At most RMAX levels of recursion are allowed, which means
2636 that at most 3**RMAX can be a factor of idim.
2637 --------------------------------------------------------------------*/
2638
2639 #undef CC3
2640 #undef SS3
2641 #define CC3 (-0.5) /* cos(2*Pi/3) */
2642 #define SS3 (0.8660254038) /* sin(2*Pi/3) */
2643
fft_3dec_OMP(int mode,int idim,complex * xc)2644 static void fft_3dec_OMP( int mode , int idim , complex *xc )
2645 {
2646 DECLARE_ithr ;
2647 int rec , mold ;
2648
2649 int N=idim , M=idim/3 , M2=2*M ;
2650 complex *cs=NULL , *aa , *bb , *cc ;
2651 register int k , i ;
2652 register float aar,aai, tr,ti, bbr,bbi, ccr,cci,
2653 t1,t2,t4,t5,t6,t8 ;
2654
2655 /* ININFO_message("fft_3dec_OMP(%d)",idim) ; */
2656
2657 /*-- initialize cosine/sine table and memory space --*/
2658
2659 rec = CFO_recX[ithr] ; /* recursion level (0 .. RMAX-1) */
2660 if( rec >= RMAX ){
2661 fprintf(stderr,"\n*** fft_3dec_OMP: too many recursions at idim=%d\n",idim); EXIT(1);
2662 }
2663
2664 CFO_assign(cs,M2) ;
2665 CFO_assign(aa,M) ;
2666 CFO_assign(bb,M) ;
2667 CFO_assign(cc,M) ;
2668
2669 mold = CFO_rmold3[rec][ithr] ; /* what was done before at this recursion level */
2670
2671 { if( M != mold ){
2672 double th = (2.0*PI/N) ;
2673 cs[0].r = 1.0 ; cs[0].i = 0.0 ;
2674 for( k=1 ; k < M2 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); }
2675 CFO_rmold3[rec][ithr] = M ;
2676 }}
2677
2678 /*-- load subarrays, and FFT each one --*/
2679
2680 for( i=k=0; k < M; k++ ){ aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; }
2681
2682 CFO_recX[ithr]++ ;
2683 switch( M ){
2684 case 1: break ;
2685
2686 case 2: fft2 (mode,aa) ; fft2 (mode,bb) ; fft2 (mode,cc) ; break;
2687 case 4: fft4 (mode,aa) ; fft4 (mode,bb) ; fft4 (mode,cc) ; break;
2688 case 8: fft8 (mode,aa) ; fft8 (mode,bb) ; fft8 (mode,cc) ; break;
2689 case 16: fft16 (mode,aa) ; fft16 (mode,bb) ; fft16 (mode,cc) ; break;
2690 case 32: fft32 (mode,aa) ; fft32 (mode,bb) ; fft32 (mode,cc) ; break;
2691 case 64: fft64 (mode,aa) ; fft64 (mode,bb) ; fft64 (mode,cc) ; break;
2692 case 128: fft128(mode,aa) ; fft128(mode,bb) ; fft128(mode,cc) ; break;
2693 case 256: fft256(mode,aa) ; fft256(mode,bb) ; fft256(mode,cc) ; break;
2694 case 512: fft512(mode,aa) ; fft512(mode,bb) ; fft512(mode,cc) ; break;
2695
2696 case 1024: fft1024(mode,aa) ; fft1024(mode,bb) ; fft1024(mode,cc) ; break;
2697 #if 1
2698 case 2048: fft2048(mode,aa) ; fft2048(mode,bb) ; fft2048(mode,cc) ; break;
2699 case 4096: fft4096(mode,aa) ; fft4096(mode,bb) ; fft4096(mode,cc) ; break;
2700 case 8192: fft8192(mode,aa) ; fft8192(mode,bb) ; fft8192(mode,cc) ; break;
2701
2702 case 16384: fft16384(mode,aa); fft16384(mode,bb); fft16384(mode,cc); break;
2703 #endif
2704
2705 default: csfft_cox_OMP(mode,M,aa); csfft_cox_OMP(mode,M,bb); csfft_cox_OMP(mode,M,cc);
2706 break ;
2707 }
2708 CFO_recX[ithr]-- ;
2709
2710 /*-- recombination --*/
2711
2712 if( mode > 0 ){
2713 for( k=0 ; k < M ; k++ ){
2714 t1 = bb[k].r; t2 = bb[k].i;
2715 tr = cs[k].r; ti = cs[k].i;
2716 bbr = tr*t1 - ti*t2 ; bbi = tr*t2 + ti*t1 ; /* b[k]*exp(+2*Pi*k/N) */
2717
2718 t1 = cc[ k].r; t2 = cc[ k].i;
2719 tr = cs[2*k].r; ti = cs[2*k].i;
2720 ccr = tr*t1 - ti*t2 ; cci = tr*t2 + ti*t1 ; /* c[k]*exp(+4*Pi*k/N) */
2721
2722 t4 = bbr+ccr ; t1 = t4*CC3 ;
2723 t8 = bbi+cci ; t6 = t8*CC3 ;
2724 t5 = (bbr-ccr)*SS3; t2 = (bbi-cci)*SS3;
2725
2726 aar = aa[k].r; aai = aa[k].i; /* a[k] */
2727
2728 xc[k ].r = aar+t4 ; xc[k ].i = aai+t8 ;
2729 xc[k+M ].r = (aar+t1)-t2 ; xc[k+M ].i = (aai+t6)+t5 ;
2730 xc[k+M2].r = (aar+t1)+t2 ; xc[k+M2].i = (aai+t6)-t5 ;
2731 }
2732 } else {
2733 for( k=0 ; k < M ; k++ ){
2734 t1 = bb[k].r; t2 = bb[k].i;
2735 tr = cs[k].r; ti = cs[k].i;
2736 bbr = tr*t1 + ti*t2 ; bbi = tr*t2 - ti*t1 ; /* b[k]*exp(-2*Pi*k/N) */
2737
2738 t1 = cc[ k].r; t2 = cc[ k].i;
2739 tr = cs[2*k].r; ti = cs[2*k].i;
2740 ccr = tr*t1 + ti*t2 ; cci = tr*t2 - ti*t1 ; /* c[k]*exp(-4*Pi*k/N) */
2741
2742 t4 = bbr+ccr ; t1 = t4*CC3 ;
2743 t8 = bbi+cci ; t6 = t8*CC3 ;
2744 t5 = (bbr-ccr)*SS3; t2 = (bbi-cci)*SS3;
2745
2746 aar = aa[k].r; aai = aa[k].i; /* a[k] */
2747
2748 xc[k ].r = aar+t4 ; xc[k ].i = aai+t8 ;
2749 xc[k+M ].r = (aar+t1)+t2 ; xc[k+M ].i = (aai+t6)-t5 ;
2750 xc[k+M2].r = (aar+t1)-t2 ; xc[k+M2].i = (aai+t6)+t5 ;
2751 }
2752 }
2753
2754 return ;
2755 }
2756
2757 /*------------------------------------------------------------------
2758 fft_5dec: do a decimation-by-5, plus a power-of-2.
2759 At most RMAX levels of recursion are allowed, which means
2760 that at most 5**RMAX can be a factor of idim.
2761 --------------------------------------------------------------------*/
2762
2763 #undef COS72
2764 #undef SIN72
2765 #define COS72 0.30901699 /* cos(72 deg) */
2766 #define SIN72 0.95105652 /* sin(72 deg) */
2767
fft_5dec_OMP(int mode,int idim,complex * xc)2768 static void fft_5dec_OMP( int mode , int idim , complex *xc )
2769 {
2770 DECLARE_ithr ;
2771 int rec , mold ;
2772 int N=idim , M=idim/5 , M2=2*M , M3=3*M , M4=4*M ;
2773 complex *cs, *aa, *bb, *cc, *dd, *ee ;
2774 register int k , i ;
2775 register float aar,aai,bbr,bbi,ccr,cci,ddr,ddi,eer,eei ;
2776 register float akp,akm,bkp,bkm , ajp,ajm,bjp,bjm ;
2777 register float ak,bk,aj,bj ;
2778 float c72 , s72 , c2,s2 ;
2779 int ss ;
2780
2781 /* ININFO_message("fft_5dec_OMP(%d)",idim) ; */
2782
2783 /*-- initialize cosine/sine table and memory space --*/
2784
2785 rec = CFO_recX[ithr] ; /* recursion level (0 .. RMAX-1) */
2786 if( rec >= RMAX ){
2787 fprintf(stderr,"\n*** fft_5dec_OMP: too many recursions at idim=%d\n",idim); EXIT(1);
2788 }
2789
2790 CFO_assign(cs,M4) ;
2791 CFO_assign(aa,M) ;
2792 CFO_assign(bb,M) ;
2793 CFO_assign(cc,M) ;
2794 CFO_assign(dd,M) ;
2795 CFO_assign(ee,M) ;
2796
2797 mold = CFO_rmold5[rec][ithr] ; /* what was done before at this recursion level */
2798
2799 { if( M != mold ){ /* new M? */
2800 double th = (2.0*PI/N) ;
2801 cs[0].r = 1.0 ; cs[0].i = 0.0 ;
2802 for( k=1 ; k < M4 ; k++ ){ cs[k].r = cos(k*th); cs[k].i = sin(k*th); }
2803 CFO_rmold5[rec][ithr] = M ;
2804 }}
2805
2806 /*-- load subarrays, and FFT each one --*/
2807
2808 for( i=k=0; k < M; k++ ){
2809 aa[k]=xc[i++]; bb[k]=xc[i++]; cc[k]=xc[i++]; dd[k]=xc[i++]; ee[k]=xc[i++];
2810 }
2811
2812 CFO_recX[ithr]++ ;
2813 switch( M ){
2814 case 1: break ;
2815
2816 case 2: fft2 (mode,aa) ; fft2 (mode,bb) ; fft2 (mode,cc) ;
2817 fft2 (mode,dd) ; fft2 (mode,ee) ; break;
2818
2819 case 4: fft4 (mode,aa) ; fft4 (mode,bb) ; fft4 (mode,cc) ;
2820 fft4 (mode,dd) ; fft4 (mode,ee) ; break;
2821
2822 case 8: fft8 (mode,aa) ; fft8 (mode,bb) ; fft8 (mode,cc) ;
2823 fft8 (mode,dd) ; fft8 (mode,ee) ; break;
2824
2825 case 16: fft16 (mode,aa) ; fft16 (mode,bb) ; fft16 (mode,cc) ;
2826 fft16 (mode,dd) ; fft16 (mode,ee) ; break;
2827
2828 case 32: fft32 (mode,aa) ; fft32 (mode,bb) ; fft32 (mode,cc) ;
2829 fft32 (mode,dd) ; fft32 (mode,ee) ; break;
2830
2831 case 64: fft64 (mode,aa) ; fft64 (mode,bb) ; fft64 (mode,cc) ;
2832 fft64 (mode,dd) ; fft64 (mode,ee) ; break;
2833
2834 case 128: fft128(mode,aa) ; fft128(mode,bb) ; fft128(mode,cc) ;
2835 fft128(mode,dd) ; fft128(mode,ee) ; break;
2836
2837 case 256: fft256(mode,aa) ; fft256(mode,bb) ; fft256(mode,cc) ;
2838 fft256(mode,dd) ; fft256(mode,ee) ; break;
2839
2840 case 512: fft512(mode,aa) ; fft512(mode,bb) ; fft512(mode,cc) ;
2841 fft512(mode,dd) ; fft512(mode,ee) ; break;
2842
2843 case 1024: fft1024(mode,aa) ; fft1024(mode,bb) ; fft1024(mode,cc) ;
2844 fft1024(mode,dd) ; fft1024(mode,ee) ; break;
2845
2846 #if 1
2847 case 2048: fft2048(mode,aa) ; fft2048(mode,bb) ; fft2048(mode,cc) ;
2848 fft2048(mode,dd) ; fft2048(mode,ee) ; break;
2849
2850 case 4096: fft4096(mode,aa) ; fft4096(mode,bb) ; fft4096(mode,cc) ;
2851 fft4096(mode,dd) ; fft4096(mode,ee) ; break;
2852
2853 case 8192: fft8192(mode,aa) ; fft8192(mode,bb) ; fft8192(mode,cc) ;
2854 fft8192(mode,dd) ; fft8192(mode,ee) ; break;
2855
2856 case 16384: fft16384(mode,aa) ; fft16384(mode,bb) ; fft16384(mode,cc) ;
2857 fft16384(mode,dd) ; fft16384(mode,ee) ; break;
2858 #endif
2859
2860 default: csfft_cox_OMP(mode,M,aa) ;
2861 csfft_cox_OMP(mode,M,bb) ; csfft_cox_OMP(mode,M,cc) ;
2862 csfft_cox_OMP(mode,M,dd) ; csfft_cox_OMP(mode,M,ee) ; break ;
2863 }
2864 CFO_recX[ithr]-- ;
2865
2866 /*-- recombination --*/
2867
2868 if( mode > 0 ){ s72 = SIN72 ; ss = 1 ; }
2869 else { s72 = -SIN72 ; ss = -1 ; }
2870
2871 c72 = COS72 ;
2872 c2 = c72 * c72 - s72 * s72 ;
2873 s2 = 2.0 * c72 * s72 ;
2874
2875 for( k=0 ; k < M ; k++ ){
2876 ak = bb[k].r ; bk = bb[k].i ;
2877 aj = cs[k].r ; bj = cs[k].i * ss ;
2878 bbr = aj*ak - bj*bk ; bbi = aj*bk + bj*ak ; /* b[k]*exp(ss*2*Pi*k/N) */
2879
2880 ak = cc[ k].r ; bk = cc[ k].i ;
2881 aj = cs[2*k].r ; bj = cs[2*k].i * ss ;
2882 ccr = aj*ak - bj*bk ; cci = aj*bk + bj*ak ; /* c[k]*exp(ss*4*Pi*k/N) */
2883
2884 ak = dd[ k].r ; bk = dd[ k].i ;
2885 aj = cs[3*k].r ; bj = cs[3*k].i * ss ;
2886 ddr = aj*ak - bj*bk ; ddi = aj*bk + bj*ak ; /* d[k]*exp(ss*6*Pi*k/N) */
2887
2888 ak = ee[ k].r ; bk = ee[ k].i ;
2889 aj = cs[4*k].r ; bj = cs[4*k].i * ss ;
2890 eer = aj*ak - bj*bk ; eei = aj*bk + bj*ak ; /* e[k]*exp(ss*8*Pi*k/N) */
2891
2892 aar = aa[k].r ; aai = aa[k].i ; /* a[k] */
2893
2894 /* the code below is (heavily) adapted from fftn.c */
2895
2896 akp = bbr + eer ; akm = bbr - eer ;
2897 bkp = bbi + eei ; bkm = bbi - eei ;
2898 ajp = ccr + ddr ; ajm = ccr - ddr ;
2899 bjp = cci + ddi ; bjm = cci - ddi ;
2900
2901 xc[k].r = aar + akp + ajp ; xc[k].i = aai + bkp + bjp ;
2902
2903 ak = akp * c72 + ajp * c2 + aar ; bk = bkp * c72 + bjp * c2 + aai ;
2904 aj = akm * s72 + ajm * s2 ; bj = bkm * s72 + bjm * s2 ;
2905
2906 xc[k+M ].r = ak - bj ; xc[k+M ].i = bk + aj ;
2907 xc[k+M4].r = ak + bj ; xc[k+M4].i = bk - aj ;
2908
2909 ak = akp * c2 + ajp * c72 + aar ; bk = bkp * c2 + bjp * c72 + aai ;
2910 aj = akm * s2 - ajm * s72 ; bj = bkm * s2 - bjm * s72 ;
2911
2912 xc[k+M2].r = ak - bj ; xc[k+M2].i = bk + aj ;
2913 xc[k+M3].r = ak + bj ; xc[k+M3].i = bk - aj ;
2914 }
2915
2916 return ;
2917 }
2918