tg2sip/libtgvoip/webrtc_dsp/modules/third_party/fft/fft.c

943 lines
26 KiB
C
Raw Normal View History

Squashed 'libtgvoip/' changes from 6053cf5..cfd62e6 cfd62e6 Why did it change the OS X project 3a58a16 2.4.3 c4a48b3 Updated OS X project 564eada Fix #63 4f64e2e fixes 0c732e2 fixes 12e76ed better logging f015b79 Merge pull request #62 from xvitaly/big-endian a1df90f Set preferred audio session parameters on iOS 59a975b Fixes 8fd89fc Fixes, mic level testing and volume adjustment 243acfa Backported WebRTC upstream patch with Big Endian support. fed3bb7 Detect when proxy does not support UDP and persist that across calls a7546d4 Merge commit '6d03dd9ae4bf48d7344341cdd2d055ebd3a6a42e' into public 6d03dd9 version 69adf70 Use server config for APM + iOS crash fix 0b42ec8 Update iOS project f1b9e63 packet logging beeea45 I apparently still suck at C++ memory management 24fceba Update project 7f54b91 crash fix f85ce99 Save more data in data saving mode f4c4f79 Collect packet stats and accept json string for server config 78e584c New protocol version: optimized packet size 8cf9177 Fixed build on iOS 9dd089d fixed build on android 5caaaaf Updated WebRTC APM cc0cf35 fixed deadlock 02f4835 Rearranged VoIPController methods and added sections 912f73d Updated OS X project 39376df Fixed audio glitches on Windows dfe1f03 Updated project 81daf3f fix 296187a Merge pull request #58 from telegramdesktop/tdesktop 44956ac Merge pull request #57 from UnigramDev/public fb0a2b0 Fix build for Linux. d6cf1b7 Updated UWP wrapper 0f06289 Merge branch 'public' of github.com:grishka/libtgvoip into public dcfad91 Fix #54 162f447 Merge pull request #56 from telegramdesktop/tdesktop a7ee511 Merge remote-tracking branch 'origin/tdesktop' into HEAD 467b148 Removed unused files b1a0b3d 2.3 9b292fd Fix warning in Xcode 10. 8d8522a Merge pull request #53 from UnigramDev/public 646f7d6 Merge branch 'public' into public 14d782b Fixes 68acf59 Added GetSignalBarsCount and GetConnectionState to CXWrapper 761c586 Added GetStats to CXWrapper f643b02 Prevent crash if UWP WASAPI devices aren't found b2ac10e Fixed UWP project 9a1ec51 Fixed build for Windows Phone, fixed some warnings 4aea54f fix git-subtree-dir: libtgvoip git-subtree-split: cfd62e66a825348ac51f49e5d20bf8827fef7a38
2019-02-06 18:22:38 +00:00
/*
* Copyright(c)1995,97 Mark Olesen <olesen@me.QueensU.CA>
* Queen's Univ at Kingston (Canada)
*
* Permission to use, copy, modify, and distribute this software for
* any purpose without fee is hereby granted, provided that this
* entire notice is included in all copies of any software which is
* or includes a copy or modification of this software and in all
* copies of the supporting documentation for such software.
*
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR
* IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR QUEEN'S
* UNIVERSITY AT KINGSTON MAKES ANY REPRESENTATION OR WARRANTY OF ANY
* KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS
* FITNESS FOR ANY PARTICULAR PURPOSE.
*
* All of which is to say that you can do what you like with this
* source code provided you don't try to sell it as your own and you
* include an unaltered copy of this message (including the
* copyright).
*
* It is also implicitly understood that bug fixes and improvements
* should make their way back to the general Internet community so
* that everyone benefits.
*
* Changes:
* Trivial type modifications by the WebRTC authors.
*/
/*
* File:
* WebRtcIsac_Fftn.c
*
* Public:
* WebRtcIsac_Fftn / fftnf ();
*
* Private:
* WebRtcIsac_Fftradix / fftradixf ();
*
* Descript:
* multivariate complex Fourier transform, computed in place
* using mixed-radix Fast Fourier Transform algorithm.
*
* Fortran code by:
* RC Singleton, Stanford Research Institute, Sept. 1968
*
* translated by f2c (version 19950721).
*
* int WebRtcIsac_Fftn (int ndim, const int dims[], REAL Re[], REAL Im[],
* int iSign, double scaling);
*
* NDIM = the total number dimensions
* DIMS = a vector of array sizes
* if NDIM is zero then DIMS must be zero-terminated
*
* RE and IM hold the real and imaginary components of the data, and return
* the resulting real and imaginary Fourier coefficients. Multidimensional
* data *must* be allocated contiguously. There is no limit on the number
* of dimensions.
*
* ISIGN = the sign of the complex exponential (ie, forward or inverse FFT)
* the magnitude of ISIGN (normally 1) is used to determine the
* correct indexing increment (see below).
*
* SCALING = normalizing constant by which the final result is *divided*
* if SCALING == -1, normalize by total dimension of the transform
* if SCALING < -1, normalize by the square-root of the total dimension
*
* example:
* tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3]
*
* int dims[3] = {n1,n2,n3}
* WebRtcIsac_Fftn (3, dims, Re, Im, 1, scaling);
*
*-----------------------------------------------------------------------*
* int WebRtcIsac_Fftradix (REAL Re[], REAL Im[], size_t nTotal, size_t nPass,
* size_t nSpan, int iSign, size_t max_factors,
* size_t max_perm);
*
* RE, IM - see above documentation
*
* Although there is no limit on the number of dimensions, WebRtcIsac_Fftradix() must
* be called once for each dimension, but the calls may be in any order.
*
* NTOTAL = the total number of complex data values
* NPASS = the dimension of the current variable
* NSPAN/NPASS = the spacing of consecutive data values while indexing the
* current variable
* ISIGN - see above documentation
*
* example:
* tri-variate transform with Re[n1][n2][n3], Im[n1][n2][n3]
*
* WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n1, n1, 1, maxf, maxp);
* WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n2, n1*n2, 1, maxf, maxp);
* WebRtcIsac_Fftradix (Re, Im, n1*n2*n3, n3, n1*n2*n3, 1, maxf, maxp);
*
* single-variate transform,
* NTOTAL = N = NSPAN = (number of complex data values),
*
* WebRtcIsac_Fftradix (Re, Im, n, n, n, 1, maxf, maxp);
*
* The data can also be stored in a single array with alternating real and
* imaginary parts, the magnitude of ISIGN is changed to 2 to give correct
* indexing increment, and data [0] and data [1] used to pass the initial
* addresses for the sequences of real and imaginary values,
*
* example:
* REAL data [2*NTOTAL];
* WebRtcIsac_Fftradix ( &data[0], &data[1], NTOTAL, nPass, nSpan, 2, maxf, maxp);
*
* for temporary allocation:
*
* MAX_FACTORS >= the maximum prime factor of NPASS
* MAX_PERM >= the number of prime factors of NPASS. In addition,
* if the square-free portion K of NPASS has two or more prime
* factors, then MAX_PERM >= (K-1)
*
* storage in FACTOR for a maximum of 15 prime factors of NPASS. if NPASS
* has more than one square-free factor, the product of the square-free
* factors must be <= 210 array storage for maximum prime factor of 23 the
* following two constants should agree with the array dimensions.
*
*----------------------------------------------------------------------*/
#include <stdlib.h>
#include <math.h>
#include "modules/third_party/fft/fft.h"
/* double precision routine */
static int
WebRtcIsac_Fftradix (double Re[], double Im[],
size_t nTotal, size_t nPass, size_t nSpan, int isign,
int max_factors, unsigned int max_perm,
FFTstr *fftstate);
#ifndef M_PI
# define M_PI 3.14159265358979323846264338327950288
#endif
#ifndef SIN60
# define SIN60 0.86602540378443865 /* sin(60 deg) */
# define COS72 0.30901699437494742 /* cos(72 deg) */
# define SIN72 0.95105651629515357 /* sin(72 deg) */
#endif
# define REAL double
# define FFTN WebRtcIsac_Fftn
# define FFTNS "fftn"
# define FFTRADIX WebRtcIsac_Fftradix
# define FFTRADIXS "fftradix"
int WebRtcIsac_Fftns(unsigned int ndim, const int dims[],
double Re[],
double Im[],
int iSign,
double scaling,
FFTstr *fftstate)
{
size_t nSpan, nPass, nTotal;
unsigned int i;
int ret, max_factors, max_perm;
/*
* tally the number of elements in the data array
* and determine the number of dimensions
*/
nTotal = 1;
if (ndim && dims [0])
{
for (i = 0; i < ndim; i++)
{
if (dims [i] <= 0)
{
return -1;
}
nTotal *= dims [i];
}
}
else
{
ndim = 0;
for (i = 0; dims [i]; i++)
{
if (dims [i] <= 0)
{
return -1;
}
nTotal *= dims [i];
ndim++;
}
}
/* determine maximum number of factors and permuations */
#if 1
/*
* follow John Beale's example, just use the largest dimension and don't
* worry about excess allocation. May be someone else will do it?
*/
max_factors = max_perm = 1;
for (i = 0; i < ndim; i++)
{
nSpan = dims [i];
if ((int)nSpan > max_factors)
{
max_factors = (int)nSpan;
}
if ((int)nSpan > max_perm)
{
max_perm = (int)nSpan;
}
}
#else
/* use the constants used in the original Fortran code */
max_factors = 23;
max_perm = 209;
#endif
/* loop over the dimensions: */
nPass = 1;
for (i = 0; i < ndim; i++)
{
nSpan = dims [i];
nPass *= nSpan;
ret = FFTRADIX (Re, Im, nTotal, nSpan, nPass, iSign,
max_factors, max_perm, fftstate);
/* exit, clean-up already done */
if (ret)
return ret;
}
/* Divide through by the normalizing constant: */
if (scaling && scaling != 1.0)
{
if (iSign < 0) iSign = -iSign;
if (scaling < 0.0)
{
scaling = (double)nTotal;
if (scaling < -1.0)
scaling = sqrt (scaling);
}
scaling = 1.0 / scaling; /* multiply is often faster */
for (i = 0; i < nTotal; i += iSign)
{
Re [i] *= scaling;
Im [i] *= scaling;
}
}
return 0;
}
/*
* singleton's mixed radix routine
*
* could move allocation out to WebRtcIsac_Fftn(), but leave it here so that it's
* possible to make this a standalone function
*/
static int 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)
{
int ii, mfactor, kspan, ispan, inc;
int j, jc, jf, jj, k, k1, k2, k3, k4, kk, kt, nn, ns, nt;
REAL radf;
REAL c1, c2, c3, cd, aa, aj, ak, ajm, ajp, akm, akp;
REAL s1, s2, s3, sd, bb, bj, bk, bjm, bjp, bkm, bkp;
REAL *Rtmp = NULL; /* temp space for real part*/
REAL *Itmp = NULL; /* temp space for imaginary part */
REAL *Cos = NULL; /* Cosine values */
REAL *Sin = NULL; /* Sine values */
REAL s60 = SIN60; /* sin(60 deg) */
REAL c72 = COS72; /* cos(72 deg) */
REAL s72 = SIN72; /* sin(72 deg) */
REAL pi2 = M_PI; /* use PI first, 2 PI later */
fftstate->SpaceAlloced = 0;
fftstate->MaxPermAlloced = 0;
// initialize to avoid warnings
k3 = c2 = c3 = s2 = s3 = 0.0;
if (nPass < 2)
return 0;
/* allocate storage */
if (fftstate->SpaceAlloced < max_factors * sizeof (REAL))
{
#ifdef SUN_BROKEN_REALLOC
if (!fftstate->SpaceAlloced) /* first time */
{
fftstate->SpaceAlloced = max_factors * sizeof (REAL);
}
else
{
#endif
fftstate->SpaceAlloced = max_factors * sizeof (REAL);
#ifdef SUN_BROKEN_REALLOC
}
#endif
}
else
{
/* allow full use of alloc'd space */
max_factors = fftstate->SpaceAlloced / sizeof (REAL);
}
if (fftstate->MaxPermAlloced < max_perm)
{
#ifdef SUN_BROKEN_REALLOC
if (!fftstate->MaxPermAlloced) /* first time */
else
#endif
fftstate->MaxPermAlloced = max_perm;
}
else
{
/* allow full use of alloc'd space */
max_perm = fftstate->MaxPermAlloced;
}
/* assign pointers */
Rtmp = (REAL *) fftstate->Tmp0;
Itmp = (REAL *) fftstate->Tmp1;
Cos = (REAL *) fftstate->Tmp2;
Sin = (REAL *) fftstate->Tmp3;
/*
* Function Body
*/
inc = iSign;
if (iSign < 0) {
s72 = -s72;
s60 = -s60;
pi2 = -pi2;
inc = -inc; /* absolute value */
}
/* adjust for strange increments */
nt = inc * (int)nTotal;
ns = inc * (int)nSpan;
kspan = ns;
nn = nt - inc;
jc = ns / (int)nPass;
radf = pi2 * (double) jc;
pi2 *= 2.0; /* use 2 PI from here on */
ii = 0;
jf = 0;
/* determine the factors of n */
mfactor = 0;
k = (int)nPass;
while (k % 16 == 0) {
mfactor++;
fftstate->factor [mfactor - 1] = 4;
k /= 16;
}
j = 3;
jj = 9;
do {
while (k % jj == 0) {
mfactor++;
fftstate->factor [mfactor - 1] = j;
k /= jj;
}
j += 2;
jj = j * j;
} while (jj <= k);
if (k <= 4) {
kt = mfactor;
fftstate->factor [mfactor] = k;
if (k != 1)
mfactor++;
} else {
if (k - (k / 4 << 2) == 0) {
mfactor++;
fftstate->factor [mfactor - 1] = 2;
k /= 4;
}
kt = mfactor;
j = 2;
do {
if (k % j == 0) {
mfactor++;
fftstate->factor [mfactor - 1] = j;
k /= j;
}
j = ((j + 1) / 2 << 1) + 1;
} while (j <= k);
}
if (kt) {
j = kt;
do {
mfactor++;
fftstate->factor [mfactor - 1] = fftstate->factor [j - 1];
j--;
} while (j);
}
/* test that mfactors is in range */
if (mfactor > FFT_NFACTOR)
{
return -1;
}
/* compute fourier transform */
for (;;) {
sd = radf / (double) kspan;
cd = sin(sd);
cd = 2.0 * cd * cd;
sd = sin(sd + sd);
kk = 0;
ii++;
switch (fftstate->factor [ii - 1]) {
case 2:
/* transform for factor of 2 (including rotation factor) */
kspan /= 2;
k1 = kspan + 2;
do {
do {
k2 = kk + kspan;
ak = Re [k2];
bk = Im [k2];
Re [k2] = Re [kk] - ak;
Im [k2] = Im [kk] - bk;
Re [kk] += ak;
Im [kk] += bk;
kk = k2 + kspan;
} while (kk < nn);
kk -= nn;
} while (kk < jc);
if (kk >= kspan)
goto Permute_Results_Label; /* exit infinite loop */
do {
c1 = 1.0 - cd;
s1 = sd;
do {
do {
do {
k2 = kk + kspan;
ak = Re [kk] - Re [k2];
bk = Im [kk] - Im [k2];
Re [kk] += Re [k2];
Im [kk] += Im [k2];
Re [k2] = c1 * ak - s1 * bk;
Im [k2] = s1 * ak + c1 * bk;
kk = k2 + kspan;
} while (kk < (nt-1));
k2 = kk - nt;
c1 = -c1;
kk = k1 - k2;
} while (kk > k2);
ak = c1 - (cd * c1 + sd * s1);
s1 = sd * c1 - cd * s1 + s1;
c1 = 2.0 - (ak * ak + s1 * s1);
s1 *= c1;
c1 *= ak;
kk += jc;
} while (kk < k2);
k1 += inc + inc;
kk = (k1 - kspan + 1) / 2 + jc - 1;
} while (kk < (jc + jc));
break;
case 4: /* transform for factor of 4 */
ispan = kspan;
kspan /= 4;
do {
c1 = 1.0;
s1 = 0.0;
do {
do {
k1 = kk + kspan;
k2 = k1 + kspan;
k3 = k2 + kspan;
akp = Re [kk] + Re [k2];
akm = Re [kk] - Re [k2];
ajp = Re [k1] + Re [k3];
ajm = Re [k1] - Re [k3];
bkp = Im [kk] + Im [k2];
bkm = Im [kk] - Im [k2];
bjp = Im [k1] + Im [k3];
bjm = Im [k1] - Im [k3];
Re [kk] = akp + ajp;
Im [kk] = bkp + bjp;
ajp = akp - ajp;
bjp = bkp - bjp;
if (iSign < 0) {
akp = akm + bjm;
bkp = bkm - ajm;
akm -= bjm;
bkm += ajm;
} else {
akp = akm - bjm;
bkp = bkm + ajm;
akm += bjm;
bkm -= ajm;
}
/* avoid useless multiplies */
if (s1 == 0.0) {
Re [k1] = akp;
Re [k2] = ajp;
Re [k3] = akm;
Im [k1] = bkp;
Im [k2] = bjp;
Im [k3] = bkm;
} else {
Re [k1] = akp * c1 - bkp * s1;
Re [k2] = ajp * c2 - bjp * s2;
Re [k3] = akm * c3 - bkm * s3;
Im [k1] = akp * s1 + bkp * c1;
Im [k2] = ajp * s2 + bjp * c2;
Im [k3] = akm * s3 + bkm * c3;
}
kk = k3 + kspan;
} while (kk < nt);
c2 = c1 - (cd * c1 + sd * s1);
s1 = sd * c1 - cd * s1 + s1;
c1 = 2.0 - (c2 * c2 + s1 * s1);
s1 *= c1;
c1 *= c2;
/* values of c2, c3, s2, s3 that will get used next time */
c2 = c1 * c1 - s1 * s1;
s2 = 2.0 * c1 * s1;
c3 = c2 * c1 - s2 * s1;
s3 = c2 * s1 + s2 * c1;
kk = kk - nt + jc;
} while (kk < kspan);
kk = kk - kspan + inc;
} while (kk < jc);
if (kspan == jc)
goto Permute_Results_Label; /* exit infinite loop */
break;
default:
/* transform for odd factors */
#ifdef FFT_RADIX4
return -1;
break;
#else /* FFT_RADIX4 */
k = fftstate->factor [ii - 1];
ispan = kspan;
kspan /= k;
switch (k) {
case 3: /* transform for factor of 3 (optional code) */
do {
do {
k1 = kk + kspan;
k2 = k1 + kspan;
ak = Re [kk];
bk = Im [kk];
aj = Re [k1] + Re [k2];
bj = Im [k1] + Im [k2];
Re [kk] = ak + aj;
Im [kk] = bk + bj;
ak -= 0.5 * aj;
bk -= 0.5 * bj;
aj = (Re [k1] - Re [k2]) * s60;
bj = (Im [k1] - Im [k2]) * s60;
Re [k1] = ak - bj;
Re [k2] = ak + bj;
Im [k1] = bk + aj;
Im [k2] = bk - aj;
kk = k2 + kspan;
} while (kk < (nn - 1));
kk -= nn;
} while (kk < kspan);
break;
case 5: /* transform for factor of 5 (optional code) */
c2 = c72 * c72 - s72 * s72;
s2 = 2.0 * c72 * s72;
do {
do {
k1 = kk + kspan;
k2 = k1 + kspan;
k3 = k2 + kspan;
k4 = k3 + kspan;
akp = Re [k1] + Re [k4];
akm = Re [k1] - Re [k4];
bkp = Im [k1] + Im [k4];
bkm = Im [k1] - Im [k4];
ajp = Re [k2] + Re [k3];
ajm = Re [k2] - Re [k3];
bjp = Im [k2] + Im [k3];
bjm = Im [k2] - Im [k3];
aa = Re [kk];
bb = Im [kk];
Re [kk] = aa + akp + ajp;
Im [kk] = bb + bkp + bjp;
ak = akp * c72 + ajp * c2 + aa;
bk = bkp * c72 + bjp * c2 + bb;
aj = akm * s72 + ajm * s2;
bj = bkm * s72 + bjm * s2;
Re [k1] = ak - bj;
Re [k4] = ak + bj;
Im [k1] = bk + aj;
Im [k4] = bk - aj;
ak = akp * c2 + ajp * c72 + aa;
bk = bkp * c2 + bjp * c72 + bb;
aj = akm * s2 - ajm * s72;
bj = bkm * s2 - bjm * s72;
Re [k2] = ak - bj;
Re [k3] = ak + bj;
Im [k2] = bk + aj;
Im [k3] = bk - aj;
kk = k4 + kspan;
} while (kk < (nn-1));
kk -= nn;
} while (kk < kspan);
break;
default:
if (k != jf) {
jf = k;
s1 = pi2 / (double) k;
c1 = cos(s1);
s1 = sin(s1);
if (jf > max_factors){
return -1;
}
Cos [jf - 1] = 1.0;
Sin [jf - 1] = 0.0;
j = 1;
do {
Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1;
Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1;
k--;
Cos [k - 1] = Cos [j - 1];
Sin [k - 1] = -Sin [j - 1];
j++;
} while (j < k);
}
do {
do {
k1 = kk;
k2 = kk + ispan;
ak = aa = Re [kk];
bk = bb = Im [kk];
j = 1;
k1 += kspan;
do {
k2 -= kspan;
j++;
Rtmp [j - 1] = Re [k1] + Re [k2];
ak += Rtmp [j - 1];
Itmp [j - 1] = Im [k1] + Im [k2];
bk += Itmp [j - 1];
j++;
Rtmp [j - 1] = Re [k1] - Re [k2];
Itmp [j - 1] = Im [k1] - Im [k2];
k1 += kspan;
} while (k1 < k2);
Re [kk] = ak;
Im [kk] = bk;
k1 = kk;
k2 = kk + ispan;
j = 1;
do {
k1 += kspan;
k2 -= kspan;
jj = j;
ak = aa;
bk = bb;
aj = 0.0;
bj = 0.0;
k = 1;
do {
k++;
ak += Rtmp [k - 1] * Cos [jj - 1];
bk += Itmp [k - 1] * Cos [jj - 1];
k++;
aj += Rtmp [k - 1] * Sin [jj - 1];
bj += Itmp [k - 1] * Sin [jj - 1];
jj += j;
if (jj > jf) {
jj -= jf;
}
} while (k < jf);
k = jf - j;
Re [k1] = ak - bj;
Im [k1] = bk + aj;
Re [k2] = ak + bj;
Im [k2] = bk - aj;
j++;
} while (j < k);
kk += ispan;
} while (kk < nn);
kk -= nn;
} while (kk < kspan);
break;
}
/* multiply by rotation factor (except for factors of 2 and 4) */
if (ii == mfactor)
goto Permute_Results_Label; /* exit infinite loop */
kk = jc;
do {
c2 = 1.0 - cd;
s1 = sd;
do {
c1 = c2;
s2 = s1;
kk += kspan;
do {
do {
ak = Re [kk];
Re [kk] = c2 * ak - s2 * Im [kk];
Im [kk] = s2 * ak + c2 * Im [kk];
kk += ispan;
} while (kk < nt);
ak = s1 * s2;
s2 = s1 * c2 + c1 * s2;
c2 = c1 * c2 - ak;
kk = kk - nt + kspan;
} while (kk < ispan);
c2 = c1 - (cd * c1 + sd * s1);
s1 += sd * c1 - cd * s1;
c1 = 2.0 - (c2 * c2 + s1 * s1);
s1 *= c1;
c2 *= c1;
kk = kk - ispan + jc;
} while (kk < kspan);
kk = kk - kspan + jc + inc;
} while (kk < (jc + jc));
break;
#endif /* FFT_RADIX4 */
}
}
/* permute the results to normal order---done in two stages */
/* permutation for square factors of n */
Permute_Results_Label:
fftstate->Perm [0] = ns;
if (kt) {
k = kt + kt + 1;
if (mfactor < k)
k--;
j = 1;
fftstate->Perm [k] = jc;
do {
fftstate->Perm [j] = fftstate->Perm [j - 1] / fftstate->factor [j - 1];
fftstate->Perm [k - 1] = fftstate->Perm [k] * fftstate->factor [j - 1];
j++;
k--;
} while (j < k);
k3 = fftstate->Perm [k];
kspan = fftstate->Perm [1];
kk = jc;
k2 = kspan;
j = 1;
if (nPass != nTotal) {
/* permutation for multivariate transform */
Permute_Multi_Label:
do {
do {
k = kk + jc;
do {
/* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
kk += inc;
k2 += inc;
} while (kk < (k-1));
kk += ns - jc;
k2 += ns - jc;
} while (kk < (nt-1));
k2 = k2 - nt + kspan;
kk = kk - nt + jc;
} while (k2 < (ns-1));
do {
do {
k2 -= fftstate->Perm [j - 1];
j++;
k2 = fftstate->Perm [j] + k2;
} while (k2 > fftstate->Perm [j - 1]);
j = 1;
do {
if (kk < (k2-1))
goto Permute_Multi_Label;
kk += jc;
k2 += kspan;
} while (k2 < (ns-1));
} while (kk < (ns-1));
} else {
/* permutation for single-variate transform (optional code) */
Permute_Single_Label:
do {
/* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
kk += inc;
k2 += kspan;
} while (k2 < (ns-1));
do {
do {
k2 -= fftstate->Perm [j - 1];
j++;
k2 = fftstate->Perm [j] + k2;
} while (k2 >= fftstate->Perm [j - 1]);
j = 1;
do {
if (kk < k2)
goto Permute_Single_Label;
kk += inc;
k2 += kspan;
} while (k2 < (ns-1));
} while (kk < (ns-1));
}
jc = k3;
}
if ((kt << 1) + 1 >= mfactor)
return 0;
ispan = fftstate->Perm [kt];
/* permutation for square-free factors of n */
j = mfactor - kt;
fftstate->factor [j] = 1;
do {
fftstate->factor [j - 1] *= fftstate->factor [j];
j--;
} while (j != kt);
kt++;
nn = fftstate->factor [kt - 1] - 1;
if (nn > (int) max_perm) {
return -1;
}
j = jj = 0;
for (;;) {
k = kt + 1;
k2 = fftstate->factor [kt - 1];
kk = fftstate->factor [k - 1];
j++;
if (j > nn)
break; /* exit infinite loop */
jj += kk;
while (jj >= k2) {
jj -= k2;
k2 = kk;
k++;
kk = fftstate->factor [k - 1];
jj += kk;
}
fftstate->Perm [j - 1] = jj;
}
/* determine the permutation cycles of length greater than 1 */
j = 0;
for (;;) {
do {
j++;
kk = fftstate->Perm [j - 1];
} while (kk < 0);
if (kk != j) {
do {
k = kk;
kk = fftstate->Perm [k - 1];
fftstate->Perm [k - 1] = -kk;
} while (kk != j);
k3 = kk;
} else {
fftstate->Perm [j - 1] = -j;
if (j == nn)
break; /* exit infinite loop */
}
}
max_factors *= inc;
/* reorder a and b, following the permutation cycles */
for (;;) {
j = k3 + 1;
nt -= ispan;
ii = nt - inc + 1;
if (nt < 0)
break; /* exit infinite loop */
do {
do {
j--;
} while (fftstate->Perm [j - 1] < 0);
jj = jc;
do {
kspan = jj;
if (jj > max_factors) {
kspan = max_factors;
}
jj -= kspan;
k = fftstate->Perm [j - 1];
kk = jc * k + ii + jj;
k1 = kk + kspan - 1;
k2 = 0;
do {
k2++;
Rtmp [k2 - 1] = Re [k1];
Itmp [k2 - 1] = Im [k1];
k1 -= inc;
} while (k1 != (kk-1));
do {
k1 = kk + kspan - 1;
k2 = k1 - jc * (k + fftstate->Perm [k - 1]);
k = -fftstate->Perm [k - 1];
do {
Re [k1] = Re [k2];
Im [k1] = Im [k2];
k1 -= inc;
k2 -= inc;
} while (k1 != (kk-1));
kk = k2 + 1;
} while (k != j);
k1 = kk + kspan - 1;
k2 = 0;
do {
k2++;
Re [k1] = Rtmp [k2 - 1];
Im [k1] = Itmp [k2 - 1];
k1 -= inc;
} while (k1 != (kk-1));
} while (jj);
} while (j != 1);
}
return 0; /* exit point here */
}
/* ---------------------- end-of-file (c source) ---------------------- */