• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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