xref: /aosp_15_r20/external/webrtc/modules/third_party/fft/fft.c (revision d9f758449e529ab9291ac668be2861e7a55c2422)

/*
 * Copyright(c)1995,97 Mark Olesen <[email protected]>
 *    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) ---------------------- */