1 /*
2 * Copyright(c)1995,97 Mark Olesen <olesen@me.QueensU.CA>
3 * Queen's Univ at Kingston (Canada)
4 *
5 * Permission to use, copy, modify, and distribute this software for
6 * any purpose without fee is hereby granted, provided that this
7 * entire notice is included in all copies of any software which is
8 * or includes a copy or modification of this software and in all
9 * copies of the supporting documentation for such software.
10 *
11 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR
12 * IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR QUEEN'S
13 * UNIVERSITY AT KINGSTON MAKES ANY REPRESENTATION OR WARRANTY OF ANY
14 * KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS
15 * FITNESS FOR ANY PARTICULAR PURPOSE.
16 *
17 * All of which is to say that you can do what you like with this
18 * source code provided you don't try to sell it as your own and you
19 * include an unaltered copy of this message (including the
20 * copyright).
21 *
22 * It is also implicitly understood that bug fixes and improvements
23 * should make their way back to the general Internet community so
24 * that everyone benefits.
25 *
26 * Changes:
27 * Trivial type modifications by the WebRTC authors.
28 */
29
30
31 /*
32 * File:
33 * WebRtcIsac_Fftn.c
34 *
35 * Public:
36 * WebRtcIsac_Fftn / fftnf ();
37 *
38 * Private:
39 * WebRtcIsac_Fftradix / fftradixf ();
40 *
41 * Descript:
42 * multivariate complex Fourier transform, computed in place
43 * using mixed-radix Fast Fourier Transform algorithm.
44 *
45 * Fortran code by:
46 * RC Singleton, Stanford Research Institute, Sept. 1968
47 *
48 * translated by f2c (version 19950721).
49 *
50 * int WebRtcIsac_Fftn (int ndim, const int dims[], REAL Re[], REAL Im[],
51 * int iSign, double scaling);
52 *
53 * NDIM = the total number dimensions
54 * DIMS = a vector of array sizes
55 * if NDIM is zero then DIMS must be zero-terminated
56 *
57 * RE and IM hold the real and imaginary components of the data, and return
58 * the resulting real and imaginary Fourier coefficients. Multidimensional
59 * data *must* be allocated contiguously. There is no limit on the number
60 * of dimensions.
61 *
62 * ISIGN = the sign of the complex exponential (ie, forward or inverse FFT)
63 * the magnitude of ISIGN (normally 1) is used to determine the
64 * correct indexing increment (see below).
65 *
66 * SCALING = normalizing constant by which the final result is *divided*
67 * if SCALING == -1, normalize by total dimension of the transform
68 * if SCALING < -1, normalize by the square-root of the total dimension
69 *
70 * example:
71 * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3]
72 *
73 * int dims[3] = {n1,n2,n3}
74 * WebRtcIsac_Fftn (3, dims, Re, Im, 1, scaling);
75 *
76 *-----------------------------------------------------------------------*
77 * int WebRtcIsac_Fftradix (REAL Re[], REAL Im[], size_t nTotal, size_t nPass,
78 * size_t nSpan, int iSign, size_t max_factors,
79 * size_t max_perm);
80 *
81 * RE, IM - see above documentation
82 *
83 * Although there is no limit on the number of dimensions, WebRtcIsac_Fftradix() must
84 * be called once for each dimension, but the calls may be in any order.
85 *
86 * NTOTAL = the total number of complex data values
87 * NPASS = the dimension of the current variable
88 * NSPAN/NPASS = the spacing of consecutive data values while indexing the
89 * current variable
90 * ISIGN - see above documentation
91 *
92 * example:
93 * tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3]
94 *
95 * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n1, n1, 1, maxf, maxp);
96 * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n2, n1*n2, 1, maxf, maxp);
97 * WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n3, n1*n2*n3, 1, maxf, maxp);
98 *
99 * single-variate transform,
100 * NTOTAL = N = NSPAN = (number of complex data values),
101 *
102 * WebRtcIsac_Fftradix (Re, Im, n, n, n, 1, maxf, maxp);
103 *
104 * The data can also be stored in a single array with alternating real and
105 * imaginary parts, the magnitude of ISIGN is changed to 2 to give correct
106 * indexing increment, and data [0] and data [1] used to pass the initial
107 * addresses for the sequences of real and imaginary values,
108 *
109 * example:
110 * REAL data [2*NTOTAL];
111 * WebRtcIsac_Fftradix ( &data[0], &data[1], NTOTAL, nPass, nSpan, 2, maxf, maxp);
112 *
113 * for temporary allocation:
114 *
115 * MAX_FACTORS >= the maximum prime factor of NPASS
116 * MAX_PERM >= the number of prime factors of NPASS. In addition,
117 * if the square-free portion K of NPASS has two or more prime
118 * factors, then MAX_PERM >= (K-1)
119 *
120 * storage in FACTOR for a maximum of 15 prime factors of NPASS. if NPASS
121 * has more than one square-free factor, the product of the square-free
122 * factors must be <= 210 array storage for maximum prime factor of 23 the
123 * following two constants should agree with the array dimensions.
124 *
125 *----------------------------------------------------------------------*/
126
127 #include <stdlib.h>
128 #include <math.h>
129
130 #include "modules/third_party/fft/fft.h"
131
132 /* double precision routine */
133 static int
134 WebRtcIsac_Fftradix (double Re[], double Im[],
135 size_t nTotal, size_t nPass, size_t nSpan, int isign,
136 int max_factors, unsigned int max_perm,
137 FFTstr *fftstate);
138
139
140
141 #ifndef M_PI
142 # define M_PI 3.14159265358979323846264338327950288
143 #endif
144
145 #ifndef SIN60
146 # define SIN60 0.86602540378443865 /* sin(60 deg) */
147 # define COS72 0.30901699437494742 /* cos(72 deg) */
148 # define SIN72 0.95105651629515357 /* sin(72 deg) */
149 #endif
150
151 # define REAL double
152 # define FFTN WebRtcIsac_Fftn
153 # define FFTNS "fftn"
154 # define FFTRADIX WebRtcIsac_Fftradix
155 # define FFTRADIXS "fftradix"
156
157
WebRtcIsac_Fftns(unsigned int ndim,const int dims[],double Re[],double Im[],int iSign,double scaling,FFTstr * fftstate)158 int WebRtcIsac_Fftns(unsigned int ndim, const int dims[],
159 double Re[],
160 double Im[],
161 int iSign,
162 double scaling,
163 FFTstr *fftstate)
164 {
165
166 size_t nSpan, nPass, nTotal;
167 unsigned int i;
168 int ret, max_factors, max_perm;
169
170 /*
171 * tally the number of elements in the data array
172 * and determine the number of dimensions
173 */
174 nTotal = 1;
175 if (ndim && dims [0])
176 {
177 for (i = 0; i < ndim; i++)
178 {
179 if (dims [i] <= 0)
180 {
181 return -1;
182 }
183 nTotal *= dims [i];
184 }
185 }
186 else
187 {
188 ndim = 0;
189 for (i = 0; dims [i]; i++)
190 {
191 if (dims [i] <= 0)
192 {
193 return -1;
194 }
195 nTotal *= dims [i];
196 ndim++;
197 }
198 }
199
200 /* determine maximum number of factors and permuations */
201 #if 1
202 /*
203 * follow John Beale's example, just use the largest dimension and don't
204 * worry about excess allocation. May be someone else will do it?
205 */
206 max_factors = max_perm = 1;
207 for (i = 0; i < ndim; i++)
208 {
209 nSpan = dims [i];
210 if ((int)nSpan > max_factors)
211 {
212 max_factors = (int)nSpan;
213 }
214 if ((int)nSpan > max_perm)
215 {
216 max_perm = (int)nSpan;
217 }
218 }
219 #else
220 /* use the constants used in the original Fortran code */
221 max_factors = 23;
222 max_perm = 209;
223 #endif
224 /* loop over the dimensions: */
225 nPass = 1;
226 for (i = 0; i < ndim; i++)
227 {
228 nSpan = dims [i];
229 nPass *= nSpan;
230 ret = FFTRADIX (Re, Im, nTotal, nSpan, nPass, iSign,
231 max_factors, max_perm, fftstate);
232 /* exit, clean-up already done */
233 if (ret)
234 return ret;
235 }
236
237 /* Divide through by the normalizing constant: */
238 if (scaling && scaling != 1.0)
239 {
240 if (iSign < 0) iSign = -iSign;
241 if (scaling < 0.0)
242 {
243 scaling = (double)nTotal;
244 if (scaling < -1.0)
245 scaling = sqrt (scaling);
246 }
247 scaling = 1.0 / scaling; /* multiply is often faster */
248 for (i = 0; i < nTotal; i += iSign)
249 {
250 Re [i] *= scaling;
251 Im [i] *= scaling;
252 }
253 }
254 return 0;
255 }
256
257 /*
258 * singleton's mixed radix routine
259 *
260 * could move allocation out to WebRtcIsac_Fftn(), but leave it here so that it's
261 * possible to make this a standalone function
262 */
263
FFTRADIX(REAL Re[],REAL Im[],size_t nTotal,size_t nPass,size_t nSpan,int iSign,int max_factors,unsigned int max_perm,FFTstr * fftstate)264 static int FFTRADIX (REAL Re[],
265 REAL Im[],
266 size_t nTotal,
267 size_t nPass,
268 size_t nSpan,
269 int iSign,
270 int max_factors,
271 unsigned int max_perm,
272 FFTstr *fftstate)
273 {
274 int ii, mfactor, kspan, ispan, inc;
275 int j, jc, jf, jj, k, k1, k2, k3, k4, kk, kt, nn, ns, nt;
276
277
278 REAL radf;
279 REAL c1, c2, c3, cd, aa, aj, ak, ajm, ajp, akm, akp;
280 REAL s1, s2, s3, sd, bb, bj, bk, bjm, bjp, bkm, bkp;
281
282 REAL *Rtmp = NULL; /* temp space for real part*/
283 REAL *Itmp = NULL; /* temp space for imaginary part */
284 REAL *Cos = NULL; /* Cosine values */
285 REAL *Sin = NULL; /* Sine values */
286
287 REAL s60 = SIN60; /* sin(60 deg) */
288 REAL c72 = COS72; /* cos(72 deg) */
289 REAL s72 = SIN72; /* sin(72 deg) */
290 REAL pi2 = M_PI; /* use PI first, 2 PI later */
291
292
293 fftstate->SpaceAlloced = 0;
294 fftstate->MaxPermAlloced = 0;
295
296
297 // initialize to avoid warnings
298 k3 = c2 = c3 = s2 = s3 = 0.0;
299
300 if (nPass < 2)
301 return 0;
302
303 /* allocate storage */
304 if (fftstate->SpaceAlloced < max_factors * sizeof (REAL))
305 {
306 #ifdef SUN_BROKEN_REALLOC
307 if (!fftstate->SpaceAlloced) /* first time */
308 {
309 fftstate->SpaceAlloced = max_factors * sizeof (REAL);
310 }
311 else
312 {
313 #endif
314 fftstate->SpaceAlloced = max_factors * sizeof (REAL);
315 #ifdef SUN_BROKEN_REALLOC
316 }
317 #endif
318 }
319 else
320 {
321 /* allow full use of alloc'd space */
322 max_factors = fftstate->SpaceAlloced / sizeof (REAL);
323 }
324 if (fftstate->MaxPermAlloced < max_perm)
325 {
326 #ifdef SUN_BROKEN_REALLOC
327 if (!fftstate->MaxPermAlloced) /* first time */
328 else
329 #endif
330 fftstate->MaxPermAlloced = max_perm;
331 }
332 else
333 {
334 /* allow full use of alloc'd space */
335 max_perm = fftstate->MaxPermAlloced;
336 }
337
338 /* assign pointers */
339 Rtmp = (REAL *) fftstate->Tmp0;
340 Itmp = (REAL *) fftstate->Tmp1;
341 Cos = (REAL *) fftstate->Tmp2;
342 Sin = (REAL *) fftstate->Tmp3;
343
344 /*
345 * Function Body
346 */
347 inc = iSign;
348 if (iSign < 0) {
349 s72 = -s72;
350 s60 = -s60;
351 pi2 = -pi2;
352 inc = -inc; /* absolute value */
353 }
354
355 /* adjust for strange increments */
356 nt = inc * (int)nTotal;
357 ns = inc * (int)nSpan;
358 kspan = ns;
359
360 nn = nt - inc;
361 jc = ns / (int)nPass;
362 radf = pi2 * (double) jc;
363 pi2 *= 2.0; /* use 2 PI from here on */
364
365 ii = 0;
366 jf = 0;
367 /* determine the factors of n */
368 mfactor = 0;
369 k = (int)nPass;
370 while (k % 16 == 0) {
371 mfactor++;
372 fftstate->factor [mfactor - 1] = 4;
373 k /= 16;
374 }
375 j = 3;
376 jj = 9;
377 do {
378 while (k % jj == 0) {
379 mfactor++;
380 fftstate->factor [mfactor - 1] = j;
381 k /= jj;
382 }
383 j += 2;
384 jj = j * j;
385 } while (jj <= k);
386 if (k <= 4) {
387 kt = mfactor;
388 fftstate->factor [mfactor] = k;
389 if (k != 1)
390 mfactor++;
391 } else {
392 if (k - (k / 4 << 2) == 0) {
393 mfactor++;
394 fftstate->factor [mfactor - 1] = 2;
395 k /= 4;
396 }
397 kt = mfactor;
398 j = 2;
399 do {
400 if (k % j == 0) {
401 mfactor++;
402 fftstate->factor [mfactor - 1] = j;
403 k /= j;
404 }
405 j = ((j + 1) / 2 << 1) + 1;
406 } while (j <= k);
407 }
408 if (kt) {
409 j = kt;
410 do {
411 mfactor++;
412 fftstate->factor [mfactor - 1] = fftstate->factor [j - 1];
413 j--;
414 } while (j);
415 }
416
417 /* test that mfactors is in range */
418 if (mfactor > FFT_NFACTOR)
419 {
420 return -1;
421 }
422
423 /* compute fourier transform */
424 for (;;) {
425 sd = radf / (double) kspan;
426 cd = sin(sd);
427 cd = 2.0 * cd * cd;
428 sd = sin(sd + sd);
429 kk = 0;
430 ii++;
431
432 switch (fftstate->factor [ii - 1]) {
433 case 2:
434 /* transform for factor of 2 (including rotation factor) */
435 kspan /= 2;
436 k1 = kspan + 2;
437 do {
438 do {
439 k2 = kk + kspan;
440 ak = Re [k2];
441 bk = Im [k2];
442 Re [k2] = Re [kk] - ak;
443 Im [k2] = Im [kk] - bk;
444 Re [kk] += ak;
445 Im [kk] += bk;
446 kk = k2 + kspan;
447 } while (kk < nn);
448 kk -= nn;
449 } while (kk < jc);
450 if (kk >= kspan)
451 goto Permute_Results_Label; /* exit infinite loop */
452 do {
453 c1 = 1.0 - cd;
454 s1 = sd;
455 do {
456 do {
457 do {
458 k2 = kk + kspan;
459 ak = Re [kk] - Re [k2];
460 bk = Im [kk] - Im [k2];
461 Re [kk] += Re [k2];
462 Im [kk] += Im [k2];
463 Re [k2] = c1 * ak - s1 * bk;
464 Im [k2] = s1 * ak + c1 * bk;
465 kk = k2 + kspan;
466 } while (kk < (nt-1));
467 k2 = kk - nt;
468 c1 = -c1;
469 kk = k1 - k2;
470 } while (kk > k2);
471 ak = c1 - (cd * c1 + sd * s1);
472 s1 = sd * c1 - cd * s1 + s1;
473 c1 = 2.0 - (ak * ak + s1 * s1);
474 s1 *= c1;
475 c1 *= ak;
476 kk += jc;
477 } while (kk < k2);
478 k1 += inc + inc;
479 kk = (k1 - kspan + 1) / 2 + jc - 1;
480 } while (kk < (jc + jc));
481 break;
482
483 case 4: /* transform for factor of 4 */
484 ispan = kspan;
485 kspan /= 4;
486
487 do {
488 c1 = 1.0;
489 s1 = 0.0;
490 do {
491 do {
492 k1 = kk + kspan;
493 k2 = k1 + kspan;
494 k3 = k2 + kspan;
495 akp = Re [kk] + Re [k2];
496 akm = Re [kk] - Re [k2];
497 ajp = Re [k1] + Re [k3];
498 ajm = Re [k1] - Re [k3];
499 bkp = Im [kk] + Im [k2];
500 bkm = Im [kk] - Im [k2];
501 bjp = Im [k1] + Im [k3];
502 bjm = Im [k1] - Im [k3];
503 Re [kk] = akp + ajp;
504 Im [kk] = bkp + bjp;
505 ajp = akp - ajp;
506 bjp = bkp - bjp;
507 if (iSign < 0) {
508 akp = akm + bjm;
509 bkp = bkm - ajm;
510 akm -= bjm;
511 bkm += ajm;
512 } else {
513 akp = akm - bjm;
514 bkp = bkm + ajm;
515 akm += bjm;
516 bkm -= ajm;
517 }
518 /* avoid useless multiplies */
519 if (s1 == 0.0) {
520 Re [k1] = akp;
521 Re [k2] = ajp;
522 Re [k3] = akm;
523 Im [k1] = bkp;
524 Im [k2] = bjp;
525 Im [k3] = bkm;
526 } else {
527 Re [k1] = akp * c1 - bkp * s1;
528 Re [k2] = ajp * c2 - bjp * s2;
529 Re [k3] = akm * c3 - bkm * s3;
530 Im [k1] = akp * s1 + bkp * c1;
531 Im [k2] = ajp * s2 + bjp * c2;
532 Im [k3] = akm * s3 + bkm * c3;
533 }
534 kk = k3 + kspan;
535 } while (kk < nt);
536
537 c2 = c1 - (cd * c1 + sd * s1);
538 s1 = sd * c1 - cd * s1 + s1;
539 c1 = 2.0 - (c2 * c2 + s1 * s1);
540 s1 *= c1;
541 c1 *= c2;
542 /* values of c2, c3, s2, s3 that will get used next time */
543 c2 = c1 * c1 - s1 * s1;
544 s2 = 2.0 * c1 * s1;
545 c3 = c2 * c1 - s2 * s1;
546 s3 = c2 * s1 + s2 * c1;
547 kk = kk - nt + jc;
548 } while (kk < kspan);
549 kk = kk - kspan + inc;
550 } while (kk < jc);
551 if (kspan == jc)
552 goto Permute_Results_Label; /* exit infinite loop */
553 break;
554
555 default:
556 /* transform for odd factors */
557 #ifdef FFT_RADIX4
558 return -1;
559 break;
560 #else /* FFT_RADIX4 */
561 k = fftstate->factor [ii - 1];
562 ispan = kspan;
563 kspan /= k;
564
565 switch (k) {
566 case 3: /* transform for factor of 3 (optional code) */
567 do {
568 do {
569 k1 = kk + kspan;
570 k2 = k1 + kspan;
571 ak = Re [kk];
572 bk = Im [kk];
573 aj = Re [k1] + Re [k2];
574 bj = Im [k1] + Im [k2];
575 Re [kk] = ak + aj;
576 Im [kk] = bk + bj;
577 ak -= 0.5 * aj;
578 bk -= 0.5 * bj;
579 aj = (Re [k1] - Re [k2]) * s60;
580 bj = (Im [k1] - Im [k2]) * s60;
581 Re [k1] = ak - bj;
582 Re [k2] = ak + bj;
583 Im [k1] = bk + aj;
584 Im [k2] = bk - aj;
585 kk = k2 + kspan;
586 } while (kk < (nn - 1));
587 kk -= nn;
588 } while (kk < kspan);
589 break;
590
591 case 5: /* transform for factor of 5 (optional code) */
592 c2 = c72 * c72 - s72 * s72;
593 s2 = 2.0 * c72 * s72;
594 do {
595 do {
596 k1 = kk + kspan;
597 k2 = k1 + kspan;
598 k3 = k2 + kspan;
599 k4 = k3 + kspan;
600 akp = Re [k1] + Re [k4];
601 akm = Re [k1] - Re [k4];
602 bkp = Im [k1] + Im [k4];
603 bkm = Im [k1] - Im [k4];
604 ajp = Re [k2] + Re [k3];
605 ajm = Re [k2] - Re [k3];
606 bjp = Im [k2] + Im [k3];
607 bjm = Im [k2] - Im [k3];
608 aa = Re [kk];
609 bb = Im [kk];
610 Re [kk] = aa + akp + ajp;
611 Im [kk] = bb + bkp + bjp;
612 ak = akp * c72 + ajp * c2 + aa;
613 bk = bkp * c72 + bjp * c2 + bb;
614 aj = akm * s72 + ajm * s2;
615 bj = bkm * s72 + bjm * s2;
616 Re [k1] = ak - bj;
617 Re [k4] = ak + bj;
618 Im [k1] = bk + aj;
619 Im [k4] = bk - aj;
620 ak = akp * c2 + ajp * c72 + aa;
621 bk = bkp * c2 + bjp * c72 + bb;
622 aj = akm * s2 - ajm * s72;
623 bj = bkm * s2 - bjm * s72;
624 Re [k2] = ak - bj;
625 Re [k3] = ak + bj;
626 Im [k2] = bk + aj;
627 Im [k3] = bk - aj;
628 kk = k4 + kspan;
629 } while (kk < (nn-1));
630 kk -= nn;
631 } while (kk < kspan);
632 break;
633
634 default:
635 if (k != jf) {
636 jf = k;
637 s1 = pi2 / (double) k;
638 c1 = cos(s1);
639 s1 = sin(s1);
640 if (jf > max_factors){
641 return -1;
642 }
643 Cos [jf - 1] = 1.0;
644 Sin [jf - 1] = 0.0;
645 j = 1;
646 do {
647 Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1;
648 Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1;
649 k--;
650 Cos [k - 1] = Cos [j - 1];
651 Sin [k - 1] = -Sin [j - 1];
652 j++;
653 } while (j < k);
654 }
655 do {
656 do {
657 k1 = kk;
658 k2 = kk + ispan;
659 ak = aa = Re [kk];
660 bk = bb = Im [kk];
661 j = 1;
662 k1 += kspan;
663 do {
664 k2 -= kspan;
665 j++;
666 Rtmp [j - 1] = Re [k1] + Re [k2];
667 ak += Rtmp [j - 1];
668 Itmp [j - 1] = Im [k1] + Im [k2];
669 bk += Itmp [j - 1];
670 j++;
671 Rtmp [j - 1] = Re [k1] - Re [k2];
672 Itmp [j - 1] = Im [k1] - Im [k2];
673 k1 += kspan;
674 } while (k1 < k2);
675 Re [kk] = ak;
676 Im [kk] = bk;
677 k1 = kk;
678 k2 = kk + ispan;
679 j = 1;
680 do {
681 k1 += kspan;
682 k2 -= kspan;
683 jj = j;
684 ak = aa;
685 bk = bb;
686 aj = 0.0;
687 bj = 0.0;
688 k = 1;
689 do {
690 k++;
691 ak += Rtmp [k - 1] * Cos [jj - 1];
692 bk += Itmp [k - 1] * Cos [jj - 1];
693 k++;
694 aj += Rtmp [k - 1] * Sin [jj - 1];
695 bj += Itmp [k - 1] * Sin [jj - 1];
696 jj += j;
697 if (jj > jf) {
698 jj -= jf;
699 }
700 } while (k < jf);
701 k = jf - j;
702 Re [k1] = ak - bj;
703 Im [k1] = bk + aj;
704 Re [k2] = ak + bj;
705 Im [k2] = bk - aj;
706 j++;
707 } while (j < k);
708 kk += ispan;
709 } while (kk < nn);
710 kk -= nn;
711 } while (kk < kspan);
712 break;
713 }
714
715 /* multiply by rotation factor (except for factors of 2 and 4) */
716 if (ii == mfactor)
717 goto Permute_Results_Label; /* exit infinite loop */
718 kk = jc;
719 do {
720 c2 = 1.0 - cd;
721 s1 = sd;
722 do {
723 c1 = c2;
724 s2 = s1;
725 kk += kspan;
726 do {
727 do {
728 ak = Re [kk];
729 Re [kk] = c2 * ak - s2 * Im [kk];
730 Im [kk] = s2 * ak + c2 * Im [kk];
731 kk += ispan;
732 } while (kk < nt);
733 ak = s1 * s2;
734 s2 = s1 * c2 + c1 * s2;
735 c2 = c1 * c2 - ak;
736 kk = kk - nt + kspan;
737 } while (kk < ispan);
738 c2 = c1 - (cd * c1 + sd * s1);
739 s1 += sd * c1 - cd * s1;
740 c1 = 2.0 - (c2 * c2 + s1 * s1);
741 s1 *= c1;
742 c2 *= c1;
743 kk = kk - ispan + jc;
744 } while (kk < kspan);
745 kk = kk - kspan + jc + inc;
746 } while (kk < (jc + jc));
747 break;
748 #endif /* FFT_RADIX4 */
749 }
750 }
751
752 /* permute the results to normal order---done in two stages */
753 /* permutation for square factors of n */
754 Permute_Results_Label:
755 fftstate->Perm [0] = ns;
756 if (kt) {
757 k = kt + kt + 1;
758 if (mfactor < k)
759 k--;
760 j = 1;
761 fftstate->Perm [k] = jc;
762 do {
763 fftstate->Perm [j] = fftstate->Perm [j - 1] / fftstate->factor [j - 1];
764 fftstate->Perm [k - 1] = fftstate->Perm [k] * fftstate->factor [j - 1];
765 j++;
766 k--;
767 } while (j < k);
768 k3 = fftstate->Perm [k];
769 kspan = fftstate->Perm [1];
770 kk = jc;
771 k2 = kspan;
772 j = 1;
773 if (nPass != nTotal) {
774 /* permutation for multivariate transform */
775 Permute_Multi_Label:
776 do {
777 do {
778 k = kk + jc;
779 do {
780 /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
781 ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
782 bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
783 kk += inc;
784 k2 += inc;
785 } while (kk < (k-1));
786 kk += ns - jc;
787 k2 += ns - jc;
788 } while (kk < (nt-1));
789 k2 = k2 - nt + kspan;
790 kk = kk - nt + jc;
791 } while (k2 < (ns-1));
792 do {
793 do {
794 k2 -= fftstate->Perm [j - 1];
795 j++;
796 k2 = fftstate->Perm [j] + k2;
797 } while (k2 > fftstate->Perm [j - 1]);
798 j = 1;
799 do {
800 if (kk < (k2-1))
801 goto Permute_Multi_Label;
802 kk += jc;
803 k2 += kspan;
804 } while (k2 < (ns-1));
805 } while (kk < (ns-1));
806 } else {
807 /* permutation for single-variate transform (optional code) */
808 Permute_Single_Label:
809 do {
810 /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
811 ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
812 bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
813 kk += inc;
814 k2 += kspan;
815 } while (k2 < (ns-1));
816 do {
817 do {
818 k2 -= fftstate->Perm [j - 1];
819 j++;
820 k2 = fftstate->Perm [j] + k2;
821 } while (k2 >= fftstate->Perm [j - 1]);
822 j = 1;
823 do {
824 if (kk < k2)
825 goto Permute_Single_Label;
826 kk += inc;
827 k2 += kspan;
828 } while (k2 < (ns-1));
829 } while (kk < (ns-1));
830 }
831 jc = k3;
832 }
833
834 if ((kt << 1) + 1 >= mfactor)
835 return 0;
836 ispan = fftstate->Perm [kt];
837 /* permutation for square-free factors of n */
838 j = mfactor - kt;
839 fftstate->factor [j] = 1;
840 do {
841 fftstate->factor [j - 1] *= fftstate->factor [j];
842 j--;
843 } while (j != kt);
844 kt++;
845 nn = fftstate->factor [kt - 1] - 1;
846 if (nn > (int) max_perm) {
847 return -1;
848 }
849 j = jj = 0;
850 for (;;) {
851 k = kt + 1;
852 k2 = fftstate->factor [kt - 1];
853 kk = fftstate->factor [k - 1];
854 j++;
855 if (j > nn)
856 break; /* exit infinite loop */
857 jj += kk;
858 while (jj >= k2) {
859 jj -= k2;
860 k2 = kk;
861 k++;
862 kk = fftstate->factor [k - 1];
863 jj += kk;
864 }
865 fftstate->Perm [j - 1] = jj;
866 }
867 /* determine the permutation cycles of length greater than 1 */
868 j = 0;
869 for (;;) {
870 do {
871 j++;
872 kk = fftstate->Perm [j - 1];
873 } while (kk < 0);
874 if (kk != j) {
875 do {
876 k = kk;
877 kk = fftstate->Perm [k - 1];
878 fftstate->Perm [k - 1] = -kk;
879 } while (kk != j);
880 k3 = kk;
881 } else {
882 fftstate->Perm [j - 1] = -j;
883 if (j == nn)
884 break; /* exit infinite loop */
885 }
886 }
887 max_factors *= inc;
888 /* reorder a and b, following the permutation cycles */
889 for (;;) {
890 j = k3 + 1;
891 nt -= ispan;
892 ii = nt - inc + 1;
893 if (nt < 0)
894 break; /* exit infinite loop */
895 do {
896 do {
897 j--;
898 } while (fftstate->Perm [j - 1] < 0);
899 jj = jc;
900 do {
901 kspan = jj;
902 if (jj > max_factors) {
903 kspan = max_factors;
904 }
905 jj -= kspan;
906 k = fftstate->Perm [j - 1];
907 kk = jc * k + ii + jj;
908 k1 = kk + kspan - 1;
909 k2 = 0;
910 do {
911 k2++;
912 Rtmp [k2 - 1] = Re [k1];
913 Itmp [k2 - 1] = Im [k1];
914 k1 -= inc;
915 } while (k1 != (kk-1));
916 do {
917 k1 = kk + kspan - 1;
918 k2 = k1 - jc * (k + fftstate->Perm [k - 1]);
919 k = -fftstate->Perm [k - 1];
920 do {
921 Re [k1] = Re [k2];
922 Im [k1] = Im [k2];
923 k1 -= inc;
924 k2 -= inc;
925 } while (k1 != (kk-1));
926 kk = k2 + 1;
927 } while (k != j);
928 k1 = kk + kspan - 1;
929 k2 = 0;
930 do {
931 k2++;
932 Re [k1] = Rtmp [k2 - 1];
933 Im [k1] = Itmp [k2 - 1];
934 k1 -= inc;
935 } while (k1 != (kk-1));
936 } while (jj);
937 } while (j != 1);
938 }
939 return 0; /* exit point here */
940 }
941 /* ---------------------- end-of-file (c source) ---------------------- */
942
943