• 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   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