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 #include "fft.h"
127
128 #include <stdlib.h>
129 #include <math.h>
130
131
132
133 /* double precision routine */
134 static int
135 WebRtcIsac_Fftradix (double Re[], double Im[],
136 size_t nTotal, size_t nPass, size_t nSpan, int isign,
137 int max_factors, unsigned int max_perm,
138 FFTstr *fftstate);
139
140
141
142 #ifndef M_PI
143 # define M_PI 3.14159265358979323846264338327950288
144 #endif
145
146 #ifndef SIN60
147 # define SIN60 0.86602540378443865 /* sin(60 deg) */
148 # define COS72 0.30901699437494742 /* cos(72 deg) */
149 # define SIN72 0.95105651629515357 /* sin(72 deg) */
150 #endif
151
152 # define REAL double
153 # define FFTN WebRtcIsac_Fftn
154 # define FFTNS "fftn"
155 # define FFTRADIX WebRtcIsac_Fftradix
156 # define FFTRADIXS "fftradix"
157
158
WebRtcIsac_Fftns(unsigned int ndim,const int dims[],double Re[],double Im[],int iSign,double scaling,FFTstr * fftstate)159 int WebRtcIsac_Fftns(unsigned int ndim, const int dims[],
160 double Re[],
161 double Im[],
162 int iSign,
163 double scaling,
164 FFTstr *fftstate)
165 {
166
167 size_t nSpan, nPass, nTotal;
168 unsigned int i;
169 int ret, max_factors, max_perm;
170
171 /*
172 * tally the number of elements in the data array
173 * and determine the number of dimensions
174 */
175 nTotal = 1;
176 if (ndim && dims [0])
177 {
178 for (i = 0; i < ndim; i++)
179 {
180 if (dims [i] <= 0)
181 {
182 return -1;
183 }
184 nTotal *= dims [i];
185 }
186 }
187 else
188 {
189 ndim = 0;
190 for (i = 0; dims [i]; i++)
191 {
192 if (dims [i] <= 0)
193 {
194 return -1;
195 }
196 nTotal *= dims [i];
197 ndim++;
198 }
199 }
200
201 /* determine maximum number of factors and permuations */
202 #if 1
203 /*
204 * follow John Beale's example, just use the largest dimension and don't
205 * worry about excess allocation. May be someone else will do it?
206 */
207 max_factors = max_perm = 1;
208 for (i = 0; i < ndim; i++)
209 {
210 nSpan = dims [i];
211 if ((int)nSpan > max_factors)
212 {
213 max_factors = (int)nSpan;
214 }
215 if ((int)nSpan > max_perm)
216 {
217 max_perm = (int)nSpan;
218 }
219 }
220 #else
221 /* use the constants used in the original Fortran code */
222 max_factors = 23;
223 max_perm = 209;
224 #endif
225 /* loop over the dimensions: */
226 nPass = 1;
227 for (i = 0; i < ndim; i++)
228 {
229 nSpan = dims [i];
230 nPass *= nSpan;
231 ret = FFTRADIX (Re, Im, nTotal, nSpan, nPass, iSign,
232 max_factors, max_perm, fftstate);
233 /* exit, clean-up already done */
234 if (ret)
235 return ret;
236 }
237
238 /* Divide through by the normalizing constant: */
239 if (scaling && scaling != 1.0)
240 {
241 if (iSign < 0) iSign = -iSign;
242 if (scaling < 0.0)
243 {
244 scaling = (double)nTotal;
245 if (scaling < -1.0)
246 scaling = sqrt (scaling);
247 }
248 scaling = 1.0 / scaling; /* multiply is often faster */
249 for (i = 0; i < nTotal; i += iSign)
250 {
251 Re [i] *= scaling;
252 Im [i] *= scaling;
253 }
254 }
255 return 0;
256 }
257
258 /*
259 * singleton's mixed radix routine
260 *
261 * could move allocation out to WebRtcIsac_Fftn(), but leave it here so that it's
262 * possible to make this a standalone function
263 */
264
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)265 static int FFTRADIX (REAL Re[],
266 REAL Im[],
267 size_t nTotal,
268 size_t nPass,
269 size_t nSpan,
270 int iSign,
271 int max_factors,
272 unsigned int max_perm,
273 FFTstr *fftstate)
274 {
275 int ii, mfactor, kspan, ispan, inc;
276 int j, jc, jf, jj, k, k1, k2, k3, k4, kk, kt, nn, ns, nt;
277
278
279 REAL radf;
280 REAL c1, c2, c3, cd, aa, aj, ak, ajm, ajp, akm, akp;
281 REAL s1, s2, s3, sd, bb, bj, bk, bjm, bjp, bkm, bkp;
282
283 REAL *Rtmp = NULL; /* temp space for real part*/
284 REAL *Itmp = NULL; /* temp space for imaginary part */
285 REAL *Cos = NULL; /* Cosine values */
286 REAL *Sin = NULL; /* Sine values */
287
288 REAL s60 = SIN60; /* sin(60 deg) */
289 REAL c72 = COS72; /* cos(72 deg) */
290 REAL s72 = SIN72; /* sin(72 deg) */
291 REAL pi2 = M_PI; /* use PI first, 2 PI later */
292
293
294 fftstate->SpaceAlloced = 0;
295 fftstate->MaxPermAlloced = 0;
296
297
298 // initialize to avoid warnings
299 k3 = c2 = c3 = s2 = s3 = 0.0;
300
301 if (nPass < 2)
302 return 0;
303
304 /* allocate storage */
305 if (fftstate->SpaceAlloced < max_factors * sizeof (REAL))
306 {
307 #ifdef SUN_BROKEN_REALLOC
308 if (!fftstate->SpaceAlloced) /* first time */
309 {
310 fftstate->SpaceAlloced = max_factors * sizeof (REAL);
311 }
312 else
313 {
314 #endif
315 fftstate->SpaceAlloced = max_factors * sizeof (REAL);
316 #ifdef SUN_BROKEN_REALLOC
317 }
318 #endif
319 }
320 else
321 {
322 /* allow full use of alloc'd space */
323 max_factors = fftstate->SpaceAlloced / sizeof (REAL);
324 }
325 if (fftstate->MaxPermAlloced < max_perm)
326 {
327 #ifdef SUN_BROKEN_REALLOC
328 if (!fftstate->MaxPermAlloced) /* first time */
329 else
330 #endif
331 fftstate->MaxPermAlloced = max_perm;
332 }
333 else
334 {
335 /* allow full use of alloc'd space */
336 max_perm = fftstate->MaxPermAlloced;
337 }
338
339 /* assign pointers */
340 Rtmp = (REAL *) fftstate->Tmp0;
341 Itmp = (REAL *) fftstate->Tmp1;
342 Cos = (REAL *) fftstate->Tmp2;
343 Sin = (REAL *) fftstate->Tmp3;
344
345 /*
346 * Function Body
347 */
348 inc = iSign;
349 if (iSign < 0) {
350 s72 = -s72;
351 s60 = -s60;
352 pi2 = -pi2;
353 inc = -inc; /* absolute value */
354 }
355
356 /* adjust for strange increments */
357 nt = inc * (int)nTotal;
358 ns = inc * (int)nSpan;
359 kspan = ns;
360
361 nn = nt - inc;
362 jc = ns / (int)nPass;
363 radf = pi2 * (double) jc;
364 pi2 *= 2.0; /* use 2 PI from here on */
365
366 ii = 0;
367 jf = 0;
368 /* determine the factors of n */
369 mfactor = 0;
370 k = (int)nPass;
371 while (k % 16 == 0) {
372 mfactor++;
373 fftstate->factor [mfactor - 1] = 4;
374 k /= 16;
375 }
376 j = 3;
377 jj = 9;
378 do {
379 while (k % jj == 0) {
380 mfactor++;
381 fftstate->factor [mfactor - 1] = j;
382 k /= jj;
383 }
384 j += 2;
385 jj = j * j;
386 } while (jj <= k);
387 if (k <= 4) {
388 kt = mfactor;
389 fftstate->factor [mfactor] = k;
390 if (k != 1)
391 mfactor++;
392 } else {
393 if (k - (k / 4 << 2) == 0) {
394 mfactor++;
395 fftstate->factor [mfactor - 1] = 2;
396 k /= 4;
397 }
398 kt = mfactor;
399 j = 2;
400 do {
401 if (k % j == 0) {
402 mfactor++;
403 fftstate->factor [mfactor - 1] = j;
404 k /= j;
405 }
406 j = ((j + 1) / 2 << 1) + 1;
407 } while (j <= k);
408 }
409 if (kt) {
410 j = kt;
411 do {
412 mfactor++;
413 fftstate->factor [mfactor - 1] = fftstate->factor [j - 1];
414 j--;
415 } while (j);
416 }
417
418 /* test that mfactors is in range */
419 if (mfactor > NFACTOR)
420 {
421 return -1;
422 }
423
424 /* compute fourier transform */
425 for (;;) {
426 sd = radf / (double) kspan;
427 cd = sin(sd);
428 cd = 2.0 * cd * cd;
429 sd = sin(sd + sd);
430 kk = 0;
431 ii++;
432
433 switch (fftstate->factor [ii - 1]) {
434 case 2:
435 /* transform for factor of 2 (including rotation factor) */
436 kspan /= 2;
437 k1 = kspan + 2;
438 do {
439 do {
440 k2 = kk + kspan;
441 ak = Re [k2];
442 bk = Im [k2];
443 Re [k2] = Re [kk] - ak;
444 Im [k2] = Im [kk] - bk;
445 Re [kk] += ak;
446 Im [kk] += bk;
447 kk = k2 + kspan;
448 } while (kk < nn);
449 kk -= nn;
450 } while (kk < jc);
451 if (kk >= kspan)
452 goto Permute_Results_Label; /* exit infinite loop */
453 do {
454 c1 = 1.0 - cd;
455 s1 = sd;
456 do {
457 do {
458 do {
459 k2 = kk + kspan;
460 ak = Re [kk] - Re [k2];
461 bk = Im [kk] - Im [k2];
462 Re [kk] += Re [k2];
463 Im [kk] += Im [k2];
464 Re [k2] = c1 * ak - s1 * bk;
465 Im [k2] = s1 * ak + c1 * bk;
466 kk = k2 + kspan;
467 } while (kk < (nt-1));
468 k2 = kk - nt;
469 c1 = -c1;
470 kk = k1 - k2;
471 } while (kk > k2);
472 ak = c1 - (cd * c1 + sd * s1);
473 s1 = sd * c1 - cd * s1 + s1;
474 c1 = 2.0 - (ak * ak + s1 * s1);
475 s1 *= c1;
476 c1 *= ak;
477 kk += jc;
478 } while (kk < k2);
479 k1 += inc + inc;
480 kk = (k1 - kspan + 1) / 2 + jc - 1;
481 } while (kk < (jc + jc));
482 break;
483
484 case 4: /* transform for factor of 4 */
485 ispan = kspan;
486 kspan /= 4;
487
488 do {
489 c1 = 1.0;
490 s1 = 0.0;
491 do {
492 do {
493 k1 = kk + kspan;
494 k2 = k1 + kspan;
495 k3 = k2 + kspan;
496 akp = Re [kk] + Re [k2];
497 akm = Re [kk] - Re [k2];
498 ajp = Re [k1] + Re [k3];
499 ajm = Re [k1] - Re [k3];
500 bkp = Im [kk] + Im [k2];
501 bkm = Im [kk] - Im [k2];
502 bjp = Im [k1] + Im [k3];
503 bjm = Im [k1] - Im [k3];
504 Re [kk] = akp + ajp;
505 Im [kk] = bkp + bjp;
506 ajp = akp - ajp;
507 bjp = bkp - bjp;
508 if (iSign < 0) {
509 akp = akm + bjm;
510 bkp = bkm - ajm;
511 akm -= bjm;
512 bkm += ajm;
513 } else {
514 akp = akm - bjm;
515 bkp = bkm + ajm;
516 akm += bjm;
517 bkm -= ajm;
518 }
519 /* avoid useless multiplies */
520 if (s1 == 0.0) {
521 Re [k1] = akp;
522 Re [k2] = ajp;
523 Re [k3] = akm;
524 Im [k1] = bkp;
525 Im [k2] = bjp;
526 Im [k3] = bkm;
527 } else {
528 Re [k1] = akp * c1 - bkp * s1;
529 Re [k2] = ajp * c2 - bjp * s2;
530 Re [k3] = akm * c3 - bkm * s3;
531 Im [k1] = akp * s1 + bkp * c1;
532 Im [k2] = ajp * s2 + bjp * c2;
533 Im [k3] = akm * s3 + bkm * c3;
534 }
535 kk = k3 + kspan;
536 } while (kk < nt);
537
538 c2 = c1 - (cd * c1 + sd * s1);
539 s1 = sd * c1 - cd * s1 + s1;
540 c1 = 2.0 - (c2 * c2 + s1 * s1);
541 s1 *= c1;
542 c1 *= c2;
543 /* values of c2, c3, s2, s3 that will get used next time */
544 c2 = c1 * c1 - s1 * s1;
545 s2 = 2.0 * c1 * s1;
546 c3 = c2 * c1 - s2 * s1;
547 s3 = c2 * s1 + s2 * c1;
548 kk = kk - nt + jc;
549 } while (kk < kspan);
550 kk = kk - kspan + inc;
551 } while (kk < jc);
552 if (kspan == jc)
553 goto Permute_Results_Label; /* exit infinite loop */
554 break;
555
556 default:
557 /* transform for odd factors */
558 #ifdef FFT_RADIX4
559 return -1;
560 break;
561 #else /* FFT_RADIX4 */
562 k = fftstate->factor [ii - 1];
563 ispan = kspan;
564 kspan /= k;
565
566 switch (k) {
567 case 3: /* transform for factor of 3 (optional code) */
568 do {
569 do {
570 k1 = kk + kspan;
571 k2 = k1 + kspan;
572 ak = Re [kk];
573 bk = Im [kk];
574 aj = Re [k1] + Re [k2];
575 bj = Im [k1] + Im [k2];
576 Re [kk] = ak + aj;
577 Im [kk] = bk + bj;
578 ak -= 0.5 * aj;
579 bk -= 0.5 * bj;
580 aj = (Re [k1] - Re [k2]) * s60;
581 bj = (Im [k1] - Im [k2]) * s60;
582 Re [k1] = ak - bj;
583 Re [k2] = ak + bj;
584 Im [k1] = bk + aj;
585 Im [k2] = bk - aj;
586 kk = k2 + kspan;
587 } while (kk < (nn - 1));
588 kk -= nn;
589 } while (kk < kspan);
590 break;
591
592 case 5: /* transform for factor of 5 (optional code) */
593 c2 = c72 * c72 - s72 * s72;
594 s2 = 2.0 * c72 * s72;
595 do {
596 do {
597 k1 = kk + kspan;
598 k2 = k1 + kspan;
599 k3 = k2 + kspan;
600 k4 = k3 + kspan;
601 akp = Re [k1] + Re [k4];
602 akm = Re [k1] - Re [k4];
603 bkp = Im [k1] + Im [k4];
604 bkm = Im [k1] - Im [k4];
605 ajp = Re [k2] + Re [k3];
606 ajm = Re [k2] - Re [k3];
607 bjp = Im [k2] + Im [k3];
608 bjm = Im [k2] - Im [k3];
609 aa = Re [kk];
610 bb = Im [kk];
611 Re [kk] = aa + akp + ajp;
612 Im [kk] = bb + bkp + bjp;
613 ak = akp * c72 + ajp * c2 + aa;
614 bk = bkp * c72 + bjp * c2 + bb;
615 aj = akm * s72 + ajm * s2;
616 bj = bkm * s72 + bjm * s2;
617 Re [k1] = ak - bj;
618 Re [k4] = ak + bj;
619 Im [k1] = bk + aj;
620 Im [k4] = bk - aj;
621 ak = akp * c2 + ajp * c72 + aa;
622 bk = bkp * c2 + bjp * c72 + bb;
623 aj = akm * s2 - ajm * s72;
624 bj = bkm * s2 - bjm * s72;
625 Re [k2] = ak - bj;
626 Re [k3] = ak + bj;
627 Im [k2] = bk + aj;
628 Im [k3] = bk - aj;
629 kk = k4 + kspan;
630 } while (kk < (nn-1));
631 kk -= nn;
632 } while (kk < kspan);
633 break;
634
635 default:
636 if (k != jf) {
637 jf = k;
638 s1 = pi2 / (double) k;
639 c1 = cos(s1);
640 s1 = sin(s1);
641 if (jf > max_factors){
642 return -1;
643 }
644 Cos [jf - 1] = 1.0;
645 Sin [jf - 1] = 0.0;
646 j = 1;
647 do {
648 Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1;
649 Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1;
650 k--;
651 Cos [k - 1] = Cos [j - 1];
652 Sin [k - 1] = -Sin [j - 1];
653 j++;
654 } while (j < k);
655 }
656 do {
657 do {
658 k1 = kk;
659 k2 = kk + ispan;
660 ak = aa = Re [kk];
661 bk = bb = Im [kk];
662 j = 1;
663 k1 += kspan;
664 do {
665 k2 -= kspan;
666 j++;
667 Rtmp [j - 1] = Re [k1] + Re [k2];
668 ak += Rtmp [j - 1];
669 Itmp [j - 1] = Im [k1] + Im [k2];
670 bk += Itmp [j - 1];
671 j++;
672 Rtmp [j - 1] = Re [k1] - Re [k2];
673 Itmp [j - 1] = Im [k1] - Im [k2];
674 k1 += kspan;
675 } while (k1 < k2);
676 Re [kk] = ak;
677 Im [kk] = bk;
678 k1 = kk;
679 k2 = kk + ispan;
680 j = 1;
681 do {
682 k1 += kspan;
683 k2 -= kspan;
684 jj = j;
685 ak = aa;
686 bk = bb;
687 aj = 0.0;
688 bj = 0.0;
689 k = 1;
690 do {
691 k++;
692 ak += Rtmp [k - 1] * Cos [jj - 1];
693 bk += Itmp [k - 1] * Cos [jj - 1];
694 k++;
695 aj += Rtmp [k - 1] * Sin [jj - 1];
696 bj += Itmp [k - 1] * Sin [jj - 1];
697 jj += j;
698 if (jj > jf) {
699 jj -= jf;
700 }
701 } while (k < jf);
702 k = jf - j;
703 Re [k1] = ak - bj;
704 Im [k1] = bk + aj;
705 Re [k2] = ak + bj;
706 Im [k2] = bk - aj;
707 j++;
708 } while (j < k);
709 kk += ispan;
710 } while (kk < nn);
711 kk -= nn;
712 } while (kk < kspan);
713 break;
714 }
715
716 /* multiply by rotation factor (except for factors of 2 and 4) */
717 if (ii == mfactor)
718 goto Permute_Results_Label; /* exit infinite loop */
719 kk = jc;
720 do {
721 c2 = 1.0 - cd;
722 s1 = sd;
723 do {
724 c1 = c2;
725 s2 = s1;
726 kk += kspan;
727 do {
728 do {
729 ak = Re [kk];
730 Re [kk] = c2 * ak - s2 * Im [kk];
731 Im [kk] = s2 * ak + c2 * Im [kk];
732 kk += ispan;
733 } while (kk < nt);
734 ak = s1 * s2;
735 s2 = s1 * c2 + c1 * s2;
736 c2 = c1 * c2 - ak;
737 kk = kk - nt + kspan;
738 } while (kk < ispan);
739 c2 = c1 - (cd * c1 + sd * s1);
740 s1 += sd * c1 - cd * s1;
741 c1 = 2.0 - (c2 * c2 + s1 * s1);
742 s1 *= c1;
743 c2 *= c1;
744 kk = kk - ispan + jc;
745 } while (kk < kspan);
746 kk = kk - kspan + jc + inc;
747 } while (kk < (jc + jc));
748 break;
749 #endif /* FFT_RADIX4 */
750 }
751 }
752
753 /* permute the results to normal order---done in two stages */
754 /* permutation for square factors of n */
755 Permute_Results_Label:
756 fftstate->Perm [0] = ns;
757 if (kt) {
758 k = kt + kt + 1;
759 if (mfactor < k)
760 k--;
761 j = 1;
762 fftstate->Perm [k] = jc;
763 do {
764 fftstate->Perm [j] = fftstate->Perm [j - 1] / fftstate->factor [j - 1];
765 fftstate->Perm [k - 1] = fftstate->Perm [k] * fftstate->factor [j - 1];
766 j++;
767 k--;
768 } while (j < k);
769 k3 = fftstate->Perm [k];
770 kspan = fftstate->Perm [1];
771 kk = jc;
772 k2 = kspan;
773 j = 1;
774 if (nPass != nTotal) {
775 /* permutation for multivariate transform */
776 Permute_Multi_Label:
777 do {
778 do {
779 k = kk + jc;
780 do {
781 /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
782 ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
783 bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
784 kk += inc;
785 k2 += inc;
786 } while (kk < (k-1));
787 kk += ns - jc;
788 k2 += ns - jc;
789 } while (kk < (nt-1));
790 k2 = k2 - nt + kspan;
791 kk = kk - nt + jc;
792 } while (k2 < (ns-1));
793 do {
794 do {
795 k2 -= fftstate->Perm [j - 1];
796 j++;
797 k2 = fftstate->Perm [j] + k2;
798 } while (k2 > fftstate->Perm [j - 1]);
799 j = 1;
800 do {
801 if (kk < (k2-1))
802 goto Permute_Multi_Label;
803 kk += jc;
804 k2 += kspan;
805 } while (k2 < (ns-1));
806 } while (kk < (ns-1));
807 } else {
808 /* permutation for single-variate transform (optional code) */
809 Permute_Single_Label:
810 do {
811 /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
812 ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
813 bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
814 kk += inc;
815 k2 += kspan;
816 } while (k2 < (ns-1));
817 do {
818 do {
819 k2 -= fftstate->Perm [j - 1];
820 j++;
821 k2 = fftstate->Perm [j] + k2;
822 } while (k2 >= fftstate->Perm [j - 1]);
823 j = 1;
824 do {
825 if (kk < k2)
826 goto Permute_Single_Label;
827 kk += inc;
828 k2 += kspan;
829 } while (k2 < (ns-1));
830 } while (kk < (ns-1));
831 }
832 jc = k3;
833 }
834
835 if ((kt << 1) + 1 >= mfactor)
836 return 0;
837 ispan = fftstate->Perm [kt];
838 /* permutation for square-free factors of n */
839 j = mfactor - kt;
840 fftstate->factor [j] = 1;
841 do {
842 fftstate->factor [j - 1] *= fftstate->factor [j];
843 j--;
844 } while (j != kt);
845 kt++;
846 nn = fftstate->factor [kt - 1] - 1;
847 if (nn > (int) max_perm) {
848 return -1;
849 }
850 j = jj = 0;
851 for (;;) {
852 k = kt + 1;
853 k2 = fftstate->factor [kt - 1];
854 kk = fftstate->factor [k - 1];
855 j++;
856 if (j > nn)
857 break; /* exit infinite loop */
858 jj += kk;
859 while (jj >= k2) {
860 jj -= k2;
861 k2 = kk;
862 k++;
863 kk = fftstate->factor [k - 1];
864 jj += kk;
865 }
866 fftstate->Perm [j - 1] = jj;
867 }
868 /* determine the permutation cycles of length greater than 1 */
869 j = 0;
870 for (;;) {
871 do {
872 j++;
873 kk = fftstate->Perm [j - 1];
874 } while (kk < 0);
875 if (kk != j) {
876 do {
877 k = kk;
878 kk = fftstate->Perm [k - 1];
879 fftstate->Perm [k - 1] = -kk;
880 } while (kk != j);
881 k3 = kk;
882 } else {
883 fftstate->Perm [j - 1] = -j;
884 if (j == nn)
885 break; /* exit infinite loop */
886 }
887 }
888 max_factors *= inc;
889 /* reorder a and b, following the permutation cycles */
890 for (;;) {
891 j = k3 + 1;
892 nt -= ispan;
893 ii = nt - inc + 1;
894 if (nt < 0)
895 break; /* exit infinite loop */
896 do {
897 do {
898 j--;
899 } while (fftstate->Perm [j - 1] < 0);
900 jj = jc;
901 do {
902 kspan = jj;
903 if (jj > max_factors) {
904 kspan = max_factors;
905 }
906 jj -= kspan;
907 k = fftstate->Perm [j - 1];
908 kk = jc * k + ii + jj;
909 k1 = kk + kspan - 1;
910 k2 = 0;
911 do {
912 k2++;
913 Rtmp [k2 - 1] = Re [k1];
914 Itmp [k2 - 1] = Im [k1];
915 k1 -= inc;
916 } while (k1 != (kk-1));
917 do {
918 k1 = kk + kspan - 1;
919 k2 = k1 - jc * (k + fftstate->Perm [k - 1]);
920 k = -fftstate->Perm [k - 1];
921 do {
922 Re [k1] = Re [k2];
923 Im [k1] = Im [k2];
924 k1 -= inc;
925 k2 -= inc;
926 } while (k1 != (kk-1));
927 kk = k2 + 1;
928 } while (k != j);
929 k1 = kk + kspan - 1;
930 k2 = 0;
931 do {
932 k2++;
933 Re [k1] = Rtmp [k2 - 1];
934 Im [k1] = Itmp [k2 - 1];
935 k1 -= inc;
936 } while (k1 != (kk-1));
937 } while (jj);
938 } while (j != 1);
939 }
940 return 0; /* exit point here */
941 }
942 /* ---------------------- end-of-file (c source) ---------------------- */
943
944