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 if (fftstate->Tmp0 == NULL || fftstate->Tmp1 == NULL || fftstate->Tmp2 == NULL || fftstate->Tmp3 == NULL
339 || fftstate->Perm == NULL) {
340 return -1;
341 }
342
343 /* assign pointers */
344 Rtmp = (REAL *) fftstate->Tmp0;
345 Itmp = (REAL *) fftstate->Tmp1;
346 Cos = (REAL *) fftstate->Tmp2;
347 Sin = (REAL *) fftstate->Tmp3;
348
349 /*
350 * Function Body
351 */
352 inc = iSign;
353 if (iSign < 0) {
354 s72 = -s72;
355 s60 = -s60;
356 pi2 = -pi2;
357 inc = -inc; /* absolute value */
358 }
359
360 /* adjust for strange increments */
361 nt = inc * (int)nTotal;
362 ns = inc * (int)nSpan;
363 kspan = ns;
364
365 nn = nt - inc;
366 jc = ns / (int)nPass;
367 radf = pi2 * (double) jc;
368 pi2 *= 2.0; /* use 2 PI from here on */
369
370 ii = 0;
371 jf = 0;
372 /* determine the factors of n */
373 mfactor = 0;
374 k = (int)nPass;
375 while (k % 16 == 0) {
376 mfactor++;
377 fftstate->factor [mfactor - 1] = 4;
378 k /= 16;
379 }
380 j = 3;
381 jj = 9;
382 do {
383 while (k % jj == 0) {
384 mfactor++;
385 fftstate->factor [mfactor - 1] = j;
386 k /= jj;
387 }
388 j += 2;
389 jj = j * j;
390 } while (jj <= k);
391 if (k <= 4) {
392 kt = mfactor;
393 fftstate->factor [mfactor] = k;
394 if (k != 1)
395 mfactor++;
396 } else {
397 if (k - (k / 4 << 2) == 0) {
398 mfactor++;
399 fftstate->factor [mfactor - 1] = 2;
400 k /= 4;
401 }
402 kt = mfactor;
403 j = 2;
404 do {
405 if (k % j == 0) {
406 mfactor++;
407 fftstate->factor [mfactor - 1] = j;
408 k /= j;
409 }
410 j = ((j + 1) / 2 << 1) + 1;
411 } while (j <= k);
412 }
413 if (kt) {
414 j = kt;
415 do {
416 mfactor++;
417 fftstate->factor [mfactor - 1] = fftstate->factor [j - 1];
418 j--;
419 } while (j);
420 }
421
422 /* test that mfactors is in range */
423 if (mfactor > NFACTOR)
424 {
425 return -1;
426 }
427
428 /* compute fourier transform */
429 for (;;) {
430 sd = radf / (double) kspan;
431 cd = sin(sd);
432 cd = 2.0 * cd * cd;
433 sd = sin(sd + sd);
434 kk = 0;
435 ii++;
436
437 switch (fftstate->factor [ii - 1]) {
438 case 2:
439 /* transform for factor of 2 (including rotation factor) */
440 kspan /= 2;
441 k1 = kspan + 2;
442 do {
443 do {
444 k2 = kk + kspan;
445 ak = Re [k2];
446 bk = Im [k2];
447 Re [k2] = Re [kk] - ak;
448 Im [k2] = Im [kk] - bk;
449 Re [kk] += ak;
450 Im [kk] += bk;
451 kk = k2 + kspan;
452 } while (kk < nn);
453 kk -= nn;
454 } while (kk < jc);
455 if (kk >= kspan)
456 goto Permute_Results_Label; /* exit infinite loop */
457 do {
458 c1 = 1.0 - cd;
459 s1 = sd;
460 do {
461 do {
462 do {
463 k2 = kk + kspan;
464 ak = Re [kk] - Re [k2];
465 bk = Im [kk] - Im [k2];
466 Re [kk] += Re [k2];
467 Im [kk] += Im [k2];
468 Re [k2] = c1 * ak - s1 * bk;
469 Im [k2] = s1 * ak + c1 * bk;
470 kk = k2 + kspan;
471 } while (kk < (nt-1));
472 k2 = kk - nt;
473 c1 = -c1;
474 kk = k1 - k2;
475 } while (kk > k2);
476 ak = c1 - (cd * c1 + sd * s1);
477 s1 = sd * c1 - cd * s1 + s1;
478 c1 = 2.0 - (ak * ak + s1 * s1);
479 s1 *= c1;
480 c1 *= ak;
481 kk += jc;
482 } while (kk < k2);
483 k1 += inc + inc;
484 kk = (k1 - kspan + 1) / 2 + jc - 1;
485 } while (kk < (jc + jc));
486 break;
487
488 case 4: /* transform for factor of 4 */
489 ispan = kspan;
490 kspan /= 4;
491
492 do {
493 c1 = 1.0;
494 s1 = 0.0;
495 do {
496 do {
497 k1 = kk + kspan;
498 k2 = k1 + kspan;
499 k3 = k2 + kspan;
500 akp = Re [kk] + Re [k2];
501 akm = Re [kk] - Re [k2];
502 ajp = Re [k1] + Re [k3];
503 ajm = Re [k1] - Re [k3];
504 bkp = Im [kk] + Im [k2];
505 bkm = Im [kk] - Im [k2];
506 bjp = Im [k1] + Im [k3];
507 bjm = Im [k1] - Im [k3];
508 Re [kk] = akp + ajp;
509 Im [kk] = bkp + bjp;
510 ajp = akp - ajp;
511 bjp = bkp - bjp;
512 if (iSign < 0) {
513 akp = akm + bjm;
514 bkp = bkm - ajm;
515 akm -= bjm;
516 bkm += ajm;
517 } else {
518 akp = akm - bjm;
519 bkp = bkm + ajm;
520 akm += bjm;
521 bkm -= ajm;
522 }
523 /* avoid useless multiplies */
524 if (s1 == 0.0) {
525 Re [k1] = akp;
526 Re [k2] = ajp;
527 Re [k3] = akm;
528 Im [k1] = bkp;
529 Im [k2] = bjp;
530 Im [k3] = bkm;
531 } else {
532 Re [k1] = akp * c1 - bkp * s1;
533 Re [k2] = ajp * c2 - bjp * s2;
534 Re [k3] = akm * c3 - bkm * s3;
535 Im [k1] = akp * s1 + bkp * c1;
536 Im [k2] = ajp * s2 + bjp * c2;
537 Im [k3] = akm * s3 + bkm * c3;
538 }
539 kk = k3 + kspan;
540 } while (kk < nt);
541
542 c2 = c1 - (cd * c1 + sd * s1);
543 s1 = sd * c1 - cd * s1 + s1;
544 c1 = 2.0 - (c2 * c2 + s1 * s1);
545 s1 *= c1;
546 c1 *= c2;
547 /* values of c2, c3, s2, s3 that will get used next time */
548 c2 = c1 * c1 - s1 * s1;
549 s2 = 2.0 * c1 * s1;
550 c3 = c2 * c1 - s2 * s1;
551 s3 = c2 * s1 + s2 * c1;
552 kk = kk - nt + jc;
553 } while (kk < kspan);
554 kk = kk - kspan + inc;
555 } while (kk < jc);
556 if (kspan == jc)
557 goto Permute_Results_Label; /* exit infinite loop */
558 break;
559
560 default:
561 /* transform for odd factors */
562 #ifdef FFT_RADIX4
563 return -1;
564 break;
565 #else /* FFT_RADIX4 */
566 k = fftstate->factor [ii - 1];
567 ispan = kspan;
568 kspan /= k;
569
570 switch (k) {
571 case 3: /* transform for factor of 3 (optional code) */
572 do {
573 do {
574 k1 = kk + kspan;
575 k2 = k1 + kspan;
576 ak = Re [kk];
577 bk = Im [kk];
578 aj = Re [k1] + Re [k2];
579 bj = Im [k1] + Im [k2];
580 Re [kk] = ak + aj;
581 Im [kk] = bk + bj;
582 ak -= 0.5 * aj;
583 bk -= 0.5 * bj;
584 aj = (Re [k1] - Re [k2]) * s60;
585 bj = (Im [k1] - Im [k2]) * s60;
586 Re [k1] = ak - bj;
587 Re [k2] = ak + bj;
588 Im [k1] = bk + aj;
589 Im [k2] = bk - aj;
590 kk = k2 + kspan;
591 } while (kk < (nn - 1));
592 kk -= nn;
593 } while (kk < kspan);
594 break;
595
596 case 5: /* transform for factor of 5 (optional code) */
597 c2 = c72 * c72 - s72 * s72;
598 s2 = 2.0 * c72 * s72;
599 do {
600 do {
601 k1 = kk + kspan;
602 k2 = k1 + kspan;
603 k3 = k2 + kspan;
604 k4 = k3 + kspan;
605 akp = Re [k1] + Re [k4];
606 akm = Re [k1] - Re [k4];
607 bkp = Im [k1] + Im [k4];
608 bkm = Im [k1] - Im [k4];
609 ajp = Re [k2] + Re [k3];
610 ajm = Re [k2] - Re [k3];
611 bjp = Im [k2] + Im [k3];
612 bjm = Im [k2] - Im [k3];
613 aa = Re [kk];
614 bb = Im [kk];
615 Re [kk] = aa + akp + ajp;
616 Im [kk] = bb + bkp + bjp;
617 ak = akp * c72 + ajp * c2 + aa;
618 bk = bkp * c72 + bjp * c2 + bb;
619 aj = akm * s72 + ajm * s2;
620 bj = bkm * s72 + bjm * s2;
621 Re [k1] = ak - bj;
622 Re [k4] = ak + bj;
623 Im [k1] = bk + aj;
624 Im [k4] = bk - aj;
625 ak = akp * c2 + ajp * c72 + aa;
626 bk = bkp * c2 + bjp * c72 + bb;
627 aj = akm * s2 - ajm * s72;
628 bj = bkm * s2 - bjm * s72;
629 Re [k2] = ak - bj;
630 Re [k3] = ak + bj;
631 Im [k2] = bk + aj;
632 Im [k3] = bk - aj;
633 kk = k4 + kspan;
634 } while (kk < (nn-1));
635 kk -= nn;
636 } while (kk < kspan);
637 break;
638
639 default:
640 if (k != jf) {
641 jf = k;
642 s1 = pi2 / (double) k;
643 c1 = cos(s1);
644 s1 = sin(s1);
645 if (jf > max_factors){
646 return -1;
647 }
648 Cos [jf - 1] = 1.0;
649 Sin [jf - 1] = 0.0;
650 j = 1;
651 do {
652 Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1;
653 Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1;
654 k--;
655 Cos [k - 1] = Cos [j - 1];
656 Sin [k - 1] = -Sin [j - 1];
657 j++;
658 } while (j < k);
659 }
660 do {
661 do {
662 k1 = kk;
663 k2 = kk + ispan;
664 ak = aa = Re [kk];
665 bk = bb = Im [kk];
666 j = 1;
667 k1 += kspan;
668 do {
669 k2 -= kspan;
670 j++;
671 Rtmp [j - 1] = Re [k1] + Re [k2];
672 ak += Rtmp [j - 1];
673 Itmp [j - 1] = Im [k1] + Im [k2];
674 bk += Itmp [j - 1];
675 j++;
676 Rtmp [j - 1] = Re [k1] - Re [k2];
677 Itmp [j - 1] = Im [k1] - Im [k2];
678 k1 += kspan;
679 } while (k1 < k2);
680 Re [kk] = ak;
681 Im [kk] = bk;
682 k1 = kk;
683 k2 = kk + ispan;
684 j = 1;
685 do {
686 k1 += kspan;
687 k2 -= kspan;
688 jj = j;
689 ak = aa;
690 bk = bb;
691 aj = 0.0;
692 bj = 0.0;
693 k = 1;
694 do {
695 k++;
696 ak += Rtmp [k - 1] * Cos [jj - 1];
697 bk += Itmp [k - 1] * Cos [jj - 1];
698 k++;
699 aj += Rtmp [k - 1] * Sin [jj - 1];
700 bj += Itmp [k - 1] * Sin [jj - 1];
701 jj += j;
702 if (jj > jf) {
703 jj -= jf;
704 }
705 } while (k < jf);
706 k = jf - j;
707 Re [k1] = ak - bj;
708 Im [k1] = bk + aj;
709 Re [k2] = ak + bj;
710 Im [k2] = bk - aj;
711 j++;
712 } while (j < k);
713 kk += ispan;
714 } while (kk < nn);
715 kk -= nn;
716 } while (kk < kspan);
717 break;
718 }
719
720 /* multiply by rotation factor (except for factors of 2 and 4) */
721 if (ii == mfactor)
722 goto Permute_Results_Label; /* exit infinite loop */
723 kk = jc;
724 do {
725 c2 = 1.0 - cd;
726 s1 = sd;
727 do {
728 c1 = c2;
729 s2 = s1;
730 kk += kspan;
731 do {
732 do {
733 ak = Re [kk];
734 Re [kk] = c2 * ak - s2 * Im [kk];
735 Im [kk] = s2 * ak + c2 * Im [kk];
736 kk += ispan;
737 } while (kk < nt);
738 ak = s1 * s2;
739 s2 = s1 * c2 + c1 * s2;
740 c2 = c1 * c2 - ak;
741 kk = kk - nt + kspan;
742 } while (kk < ispan);
743 c2 = c1 - (cd * c1 + sd * s1);
744 s1 += sd * c1 - cd * s1;
745 c1 = 2.0 - (c2 * c2 + s1 * s1);
746 s1 *= c1;
747 c2 *= c1;
748 kk = kk - ispan + jc;
749 } while (kk < kspan);
750 kk = kk - kspan + jc + inc;
751 } while (kk < (jc + jc));
752 break;
753 #endif /* FFT_RADIX4 */
754 }
755 }
756
757 /* permute the results to normal order---done in two stages */
758 /* permutation for square factors of n */
759 Permute_Results_Label:
760 fftstate->Perm [0] = ns;
761 if (kt) {
762 k = kt + kt + 1;
763 if (mfactor < k)
764 k--;
765 j = 1;
766 fftstate->Perm [k] = jc;
767 do {
768 fftstate->Perm [j] = fftstate->Perm [j - 1] / fftstate->factor [j - 1];
769 fftstate->Perm [k - 1] = fftstate->Perm [k] * fftstate->factor [j - 1];
770 j++;
771 k--;
772 } while (j < k);
773 k3 = fftstate->Perm [k];
774 kspan = fftstate->Perm [1];
775 kk = jc;
776 k2 = kspan;
777 j = 1;
778 if (nPass != nTotal) {
779 /* permutation for multivariate transform */
780 Permute_Multi_Label:
781 do {
782 do {
783 k = kk + jc;
784 do {
785 /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
786 ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
787 bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
788 kk += inc;
789 k2 += inc;
790 } while (kk < (k-1));
791 kk += ns - jc;
792 k2 += ns - jc;
793 } while (kk < (nt-1));
794 k2 = k2 - nt + kspan;
795 kk = kk - nt + jc;
796 } while (k2 < (ns-1));
797 do {
798 do {
799 k2 -= fftstate->Perm [j - 1];
800 j++;
801 k2 = fftstate->Perm [j] + k2;
802 } while (k2 > fftstate->Perm [j - 1]);
803 j = 1;
804 do {
805 if (kk < (k2-1))
806 goto Permute_Multi_Label;
807 kk += jc;
808 k2 += kspan;
809 } while (k2 < (ns-1));
810 } while (kk < (ns-1));
811 } else {
812 /* permutation for single-variate transform (optional code) */
813 Permute_Single_Label:
814 do {
815 /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
816 ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
817 bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
818 kk += inc;
819 k2 += kspan;
820 } while (k2 < (ns-1));
821 do {
822 do {
823 k2 -= fftstate->Perm [j - 1];
824 j++;
825 k2 = fftstate->Perm [j] + k2;
826 } while (k2 >= fftstate->Perm [j - 1]);
827 j = 1;
828 do {
829 if (kk < k2)
830 goto Permute_Single_Label;
831 kk += inc;
832 k2 += kspan;
833 } while (k2 < (ns-1));
834 } while (kk < (ns-1));
835 }
836 jc = k3;
837 }
838
839 if ((kt << 1) + 1 >= mfactor)
840 return 0;
841 ispan = fftstate->Perm [kt];
842 /* permutation for square-free factors of n */
843 j = mfactor - kt;
844 fftstate->factor [j] = 1;
845 do {
846 fftstate->factor [j - 1] *= fftstate->factor [j];
847 j--;
848 } while (j != kt);
849 kt++;
850 nn = fftstate->factor [kt - 1] - 1;
851 if (nn > (int) max_perm) {
852 return -1;
853 }
854 j = jj = 0;
855 for (;;) {
856 k = kt + 1;
857 k2 = fftstate->factor [kt - 1];
858 kk = fftstate->factor [k - 1];
859 j++;
860 if (j > nn)
861 break; /* exit infinite loop */
862 jj += kk;
863 while (jj >= k2) {
864 jj -= k2;
865 k2 = kk;
866 k++;
867 kk = fftstate->factor [k - 1];
868 jj += kk;
869 }
870 fftstate->Perm [j - 1] = jj;
871 }
872 /* determine the permutation cycles of length greater than 1 */
873 j = 0;
874 for (;;) {
875 do {
876 j++;
877 kk = fftstate->Perm [j - 1];
878 } while (kk < 0);
879 if (kk != j) {
880 do {
881 k = kk;
882 kk = fftstate->Perm [k - 1];
883 fftstate->Perm [k - 1] = -kk;
884 } while (kk != j);
885 k3 = kk;
886 } else {
887 fftstate->Perm [j - 1] = -j;
888 if (j == nn)
889 break; /* exit infinite loop */
890 }
891 }
892 max_factors *= inc;
893 /* reorder a and b, following the permutation cycles */
894 for (;;) {
895 j = k3 + 1;
896 nt -= ispan;
897 ii = nt - inc + 1;
898 if (nt < 0)
899 break; /* exit infinite loop */
900 do {
901 do {
902 j--;
903 } while (fftstate->Perm [j - 1] < 0);
904 jj = jc;
905 do {
906 kspan = jj;
907 if (jj > max_factors) {
908 kspan = max_factors;
909 }
910 jj -= kspan;
911 k = fftstate->Perm [j - 1];
912 kk = jc * k + ii + jj;
913 k1 = kk + kspan - 1;
914 k2 = 0;
915 do {
916 k2++;
917 Rtmp [k2 - 1] = Re [k1];
918 Itmp [k2 - 1] = Im [k1];
919 k1 -= inc;
920 } while (k1 != (kk-1));
921 do {
922 k1 = kk + kspan - 1;
923 k2 = k1 - jc * (k + fftstate->Perm [k - 1]);
924 k = -fftstate->Perm [k - 1];
925 do {
926 Re [k1] = Re [k2];
927 Im [k1] = Im [k2];
928 k1 -= inc;
929 k2 -= inc;
930 } while (k1 != (kk-1));
931 kk = k2 + 1;
932 } while (k != j);
933 k1 = kk + kspan - 1;
934 k2 = 0;
935 do {
936 k2++;
937 Re [k1] = Rtmp [k2 - 1];
938 Im [k1] = Itmp [k2 - 1];
939 k1 -= inc;
940 } while (k1 != (kk-1));
941 } while (jj);
942 } while (j != 1);
943 }
944 return 0; /* exit point here */
945 }
946 /* ---------------------- end-of-file (c source) ---------------------- */
947
948