File Coverage

lib/PDL/FFT/fftn.c

Criterion	Covered	Total	%
statement	1207	1898	63.5
branch	240	698	34.3
condition			n/a
subroutine			n/a
pod			n/a
total	1447	2596	55.7

/*--------------------------------*-C-*---------------------------------*
 * File:
 *	fftn.c
 *
 *	multivariate complex Fourier transform, computed in place
 *	using mixed-radix Fast Fourier Transform algorithm.
 *
 *	Fortran code by:
 *	    RC Singleton, Stanford Research Institute, Sept. 1968
 *		NIST Guide to Available Math Software.
 *		Source for module FFT from package GO.
 *		Retrieved from NETLIB on Wed Jul  5 11:50:07 1995.
 *	translated by f2c (version 19950721) and with lots of cleanup
 *	to make it resemble C by:
 *	    MJ Olesen, Queen's University at Kingston, 1995-97
 */
/*{{{ Copyright: */
/*
 * Copyright(c)1995,97 Mark Olesen
 *		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.
 *----------------------------------------------------------------------*/
/*}}}*/
/*{{{ notes: */
/*
 * Public:
 *	fftn / fftnf / fftnl
 *	(these are documented in the header file)
 *
 * Private:
 *	fftradix / fftradixf
 *
 * ----------------------------------------------------------------------*
 * int fftradix (REAL Re[], REAL Im[], size_t nTotal, size_t nPass,
 *		 size_t nSpan, int iSign, size_t maxFactors,
 *		 size_t maxPerm);
 *
 * 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.
 *
 *
 * Although there is no limit on the number of dimensions, 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]
 *
 *	fftradix (Re, Im, n1*n2*n3, n1,       n1, 1, maxf, maxp);
 *	fftradix (Re, Im, n1*n2*n3, n2,    n1*n2, 1, maxf, maxp);
 *	fftradix (Re, Im, n1*n2*n3, n3, n1*n2*n3, 1, maxf, maxp);
 *
 * single-variate transform,
 *    NTOTAL = N = NSPAN = (number of complex data values),
 *
 *	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];
 *	fftradix (&data[0], &data[1], NTOTAL, nPass, nSpan, 2, maxf, maxp);
 *
 * for temporary allocation:
 *
 * MAXFACTORS	>= the maximum prime factor of NPASS
 * MAXPERM	>= the number of prime factors of NPASS.  In addition,
 *	if the square-free portion K of NPASS has two or more prime
 *	factors, then MAXPERM >= (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.
 * ----------------------------------------------------------------------*/
/*}}}*/
/*{{{ Revisions: */
/*
 * 26 July 95	John Beale
 *	- added maxf and maxp as parameters to fftradix()
 *
 * 28 July 95	Mark Olesen
 *	- cleaned-up the Fortran 66 goto spaghetti, only 3 labels remain.
 *
 *	- added fft_free() to provide some measure of control over
 *	  allocation/deallocation.
 *
 *	- added fftn() wrapper for multidimensional FFTs
 *
 *	- use -DFFT_NOFLOAT or -DFFT_NODOUBLE to avoid compiling that
 *	  precision. Note suffix `f' on the function names indicates
 *	  float precision.
 *
 *	- revised documentation
 *
 * 31 July 95	Mark Olesen
 *	- added GNU Public License
 *	- more cleanup
 *	- define SUN_BROKEN_REALLOC to use malloc() instead of realloc()
 *	  on the first pass through, apparently needed for old libc
 *	- removed #error directive in favour of some code that simply
 *	  won't compile (generate an error that way)
 *
 * 1 Aug 95	Mark Olesen
 *	- define FFT_RADIX4 to only have radix 2, radix 4 transforms
 *	- made fftradix /fftradixf () static scope, just use fftn()
 *	  instead.  If you have good ideas about fixing the factors
 *	  in fftn() please do so.
 *
 * 8 Jan 95	mj olesen
 *	- fixed typo's, including one that broke scaling for scaling by
 *	  total number of matrix elements or the square root of same
 *	- removed unnecessary casts from allocations
 *
 * 10 Dec 96	mj olesen
 *	- changes defines to compile *without* float support by default,
 *	  use -DFFT_FLOAT to enable.
 *	- shifted some variables to local scope	(better hints for optimizer)
 *	- added Michael Steffens
 *	  Fortran 90 module
 *	- made it simpler to pass dimensions for 1D FFT.
 *
 * 23 Feb 97	Mark Olesen
 *	- removed the GNU Public License (see 21 July 1995 entry),
 *	  which should make it clear why I have the right to do so.
 *	- Added copyright notice and submitted to netlib
 *	- Moved documentation for the public functions to the header
 *	  files where is will always be available.
 *
 * 11 Nov 97    Robin Williams
 *      - fixed all temporaries to double precision using REALFIX,
 *        following suggestion from Scott Wurcer
 * ----------------------------------------------------------------------*/
/*}}}*/
#ifndef _FFTN_C
#define _FFTN_C
/* we use CPP to re-include this same file for double/float cases */

#include
#include
#include
#include "fftn.h"

/*{{{ defines/constants */
#ifndef M_PI
# define M_PI	3.14159265358979323846264338327950288
#endif

#define FFT_STR(x) #x
#define FFT_CONCAT(x, y) FFT_CONCAT_(x, y)
#define FFT_CONCAT_(x, y) x ## y /* second call forces arg expansion */

#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 NFACTOR	11
#define REALFIX       	long double

/*
 * use macros to access the Real/Imaginary parts so that it's possible
 * to substitute different macros if a complex struct is used
 */
#define Re_Data(i)	Re[i]
#define Im_Data(i)	Im[i]

/* return the number of factors */
static int
factorize (int nPass, int * kt, int *factor)
{
   int nFactor = 0;
   int j, jj;

   *kt = 0;
   /* determine the factors of n */
   while ((nPass % 16) == 0)	/* factors of 4 */
     {
	factor [nFactor++] = 4;
	nPass /= 16;
     }
   j = 3; jj = 9;		/* factors of 3, 5, 7, ... */
   do {
      while ((nPass % jj) == 0)
	{
	   factor [nFactor++] = j;
	   nPass /= jj;
	}
      j += 2;
      jj = j * j;
   } while (jj <= nPass);
   if (nPass <= 4)
     {
	*kt = nFactor;
	factor [nFactor] = nPass;
	if (nPass != 1)
	  nFactor++;
     }
   else
     {
	if (nPass - (nPass / 4 << 2) == 0)
	  {
	     factor [nFactor++] = 2;
	     nPass /= 4;
	  }
	*kt = nFactor;
	j = 2;
	do {
	   if ((nPass % j) == 0)
	     {
		factor [nFactor++] = j;
		nPass /= j;
	     }
	   j = ((j + 1) / 2 << 1) + 1;
	} while (j <= nPass);
     }
   if (*kt)
     {
	j = *kt;
	do
	  factor [nFactor++] = factor [--j];
	while (j);
     }

   return nFactor;
}

/*{{{ defines for re-including long double precision */
#ifdef FFT_LDOUBLE
# undef REAL
# undef FFT_SUFFIX
# define FFT_SUFFIX l
# define REAL		long double
# include "fftn.c"			/* include this file again */
#endif
/*}}}*/

/*{{{ defines for re-including double precision */
#ifdef FFT_DOUBLE
# undef REAL
# undef FFT_SUFFIX
# define FFT_SUFFIX
# define REAL		double
# include "fftn.c"			/* include this file again */
#endif
/*}}}*/

/*{{{ defines for re-including float precision */
#ifdef FFT_FLOAT
# undef REAL
# undef FFT_SUFFIX
# define FFT_SUFFIX f
# define REAL		float
# include "fftn.c"			/* include this file again */
#endif
/*}}}*/
#else	/* _FFTN_C */

# undef FFTN
# undef FFTNS
# undef FFTRADIX
# undef FFTRADIXS
# define FFTN		FFT_CONCAT(fftn, FFT_SUFFIX) /* maybe trailing [fl] */
# define FFTNS		FFT_STR(FFTN)	/* name for error message */
# define FFTRADIX	FFT_CONCAT(fftradix, FFT_SUFFIX) /* maybe trailing [fl] */
# define FFTRADIXS	FFT_STR(FFTRADIX)	/* name for error message */

/*----------------------------------------------------------------------*/
/*
 * singleton's mixed radix routine
 */
static int
FFTRADIX (REAL Re [],
	  REAL Im [],
	  size_t nTotal,
	  size_t nPass,
	  size_t nSpan,
	  int iSign,
	  int maxFactors,
	  int maxPerm)
{
#ifndef FFTRADIX_RETURN
#define FFTRADIX_RETURN(v) \
  do { \
    if (Rtmp) free(Rtmp); \
    if (Itmp) free(Itmp); \
    if (Cos) free(Cos); \
    if (Sin) free(Sin); \
    if (Perm) free(Perm); \
    return v; \
  } while (0)
#endif
#ifndef FFTRADIX_MALLOC
#define FFTRADIX_MALLOC(t, p, n) \
  do { p = malloc(sizeof(t) * n); if (!p) FFTRADIX_RETURN(-1); } while (0)
#endif
   if (nPass < 2)
     return 0;

   int ii, nFactor, factor[NFACTOR], kspan, ispan, inc;
   int j, jc, jf, jj, k, k1, k3, kk, kt, nn, ns, nt;

   REALFIX radf;
   REALFIX c1, c2, c3, cd;
   REALFIX s1, s2, s3, sd;

   REALFIX *Rtmp=NULL, *Itmp=NULL, *Cos=NULL, *Sin=NULL;
   int *Perm=NULL;
   FFTRADIX_MALLOC(REALFIX, Rtmp, maxFactors); /* temp space for real part */
   FFTRADIX_MALLOC(REALFIX, Itmp, maxFactors); /* temp space for imag part */
   FFTRADIX_MALLOC(REALFIX, Cos, maxFactors); /* Cosine values */
   FFTRADIX_MALLOC(REALFIX, Sin, maxFactors); /* Sine values */
   FFTRADIX_MALLOC(int, Perm, maxPerm);

#ifndef FFT_RADIX4
   REALFIX s60 = SIN60;		/* sin(60 deg) */
   REALFIX s72 = SIN72;		/* sin(72 deg) */
   REALFIX c72 = COS72;		/* cos(72 deg) */
#endif
   REALFIX pi2 = M_PI;		/* use PI first, 2 PI later */

   /* gcc complains about k3 being uninitialized, but I can't find out where
    * or why ... it looks okay to me.
    *
    * initialize to make gcc happy
    */
   k3 = 0;

   /* gcc complains about c2, c3, s2,s3 being uninitialized, but they're
    * only used for the radix 4 case and only AFTER the (s1 == 0.0) pass
    * through the loop at which point they will have been calculated.
    *
    * initialize to make gcc happy
    */
   c2 = c3 = s2 = s3 = 0.0;

   /* Parameter adjustments, was fortran so fix zero-offset */
   Re--;
   Im--;

   /*
    * Function Body
    */
   inc = iSign;
   if (iSign < 0)
     {
#ifndef FFT_RADIX4
	s60 = -s60;
	s72 = -s72;
#endif
	pi2 = -pi2;
	inc = -inc;		/* absolute value */
     }

   /* adjust for strange increments */
   nt = inc * nTotal;
   ns = inc * nSpan;
   kspan = ns;

   nn = nt - inc;
   jc = ns / nPass;
   radf = pi2 * (double) jc;
   pi2 *= 2.0;			/* use 2 PI from here on */

   ii = 0;
   jf = 0;
   /* determine the factors of n */

   nFactor = factorize (nPass, &kt, factor);
   /* test that nFactors is in range */
   if (nFactor > NFACTOR)
     {
	fprintf (stderr, "Error: " FFTRADIXS "() - exceeded number of factors\n");
	goto Memory_Error;
     }

   /* compute fourier transform */
   for (;;) {
      sd = radf / (double) kspan;
      cd = sin (sd);
      cd = 2.0 * cd * cd;
      sd = sin (sd + sd);
      kk = 1;
      ii++;

      switch (factor [ii - 1]) {
       case 2:
	 /* transform for factor of 2 (including rotation factor) */
	 kspan /= 2;
	 k1 = kspan + 2;
	 do {
	    do {
	       REALFIX tmpr;
	       REALFIX tmpi;
	       int k2;

	       k2 = kk + kspan;
	       tmpr = Re_Data (k2);
	       tmpi = Im_Data (k2);
	       Re_Data (k2) = Re_Data (kk) - tmpr;
	       Im_Data (k2) = Im_Data (kk) - tmpi;
	       Re_Data (kk) += tmpr;
	       Im_Data (kk) += tmpi;
	       kk = k2 + kspan;
	    } while (kk <= nn);
	    kk -= nn;
	 } while (kk <= jc);
	 if (kk > kspan)
	   goto Permute_Results;	/* exit infinite loop */
	 do {
	    int k2;

	    c1 = 1.0 - cd;
	    s1 = sd;
	    do {
	       REALFIX tmp;
	       do {
		  do {
		     REALFIX tmpr;
		     REALFIX tmpi;

		     k2 = kk + kspan;
		     tmpr = Re_Data (kk) - Re_Data (k2);
		     tmpi = Im_Data (kk) - Im_Data (k2);
		     Re_Data (kk) += Re_Data (k2);
		     Im_Data (kk) += Im_Data (k2);
		     Re_Data (k2) = c1 * tmpr - s1 * tmpi;
		     Im_Data (k2) = s1 * tmpr + c1 * tmpi;
		     kk = k2 + kspan;
		  } while (kk < nt);
		  k2 = kk - nt;
		  c1 = -c1;
		  kk = k1 - k2;
	       } while (kk > k2);
	       tmp = c1 - (cd * c1 + sd * s1);
	       s1 = sd * c1 - cd * s1 + s1;
	       c1 = 2.0 - (tmp * tmp + s1 * s1);
	       s1 *= c1;
	       c1 *= tmp;
	       kk += jc;
	    } while (kk < k2);
	    k1 += (inc + inc);
	    kk = (k1 - kspan) / 2 + jc;
	 } 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 {
		  REALFIX ajm, ajp, akm, akp;
		  REALFIX bjm, bjp, bkm, bkp;
		  int k2;

		  k1 = kk + kspan;
		  k2 = k1 + kspan;
		  k3 = k2 + kspan;
		  akp = Re_Data (kk) + Re_Data (k2);
		  akm = Re_Data (kk) - Re_Data (k2);

		  ajp = Re_Data (k1) + Re_Data (k3);
		  ajm = Re_Data (k1) - Re_Data (k3);

		  bkp = Im_Data (kk) + Im_Data (k2);
		  bkm = Im_Data (kk) - Im_Data (k2);

		  bjp = Im_Data (k1) + Im_Data (k3);
		  bjm = Im_Data (k1) - Im_Data (k3);

		  Re_Data (kk) = akp + ajp;
		  Im_Data (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_Data (k1) = akp;
		       Re_Data (k2) = ajp;
		       Re_Data (k3) = akm;
		       Im_Data (k1) = bkp;
		       Im_Data (k2) = bjp;
		       Im_Data (k3) = bkm;
		    }
		  else
		    {
		       Re_Data (k1) = akp * c1 - bkp * s1;
		       Re_Data (k2) = ajp * c2 - bjp * s2;
		       Re_Data (k3) = akm * c3 - bkm * s3;
		       Im_Data (k1) = akp * s1 + bkp * c1;
		       Im_Data (k2) = ajp * s2 + bjp * c2;
		       Im_Data (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;	/* exit infinite loop */
	 break;

       default:
	 /* transform for odd factors */
#ifdef FFT_RADIX4
	 fprintf (stderr, "Error: " FFTRADIXS "(): compiled for radix 2/4 only\n");
	 FFTRADIX_RETURN(-1);
	 break;
#else	/* FFT_RADIX4 */
	 ispan = kspan;
	 k = factor [ii - 1];
	 kspan /= factor [ii - 1];

	 switch (factor [ii - 1]) {
	  case 3:	/* transform for factor of 3 (optional code) */
	    do {
	       do {
		  REALFIX aj, tmpr;
		  REALFIX bj, tmpi;
		  int k2;

		  k1 = kk + kspan;
		  k2 = k1 + kspan;
		  tmpr = Re_Data (kk);
		  tmpi = Im_Data (kk);
		  aj = Re_Data (k1) + Re_Data (k2);
		  bj = Im_Data (k1) + Im_Data (k2);
		  Re_Data (kk) = tmpr + aj;
		  Im_Data (kk) = tmpi + bj;
		  tmpr -= 0.5 * aj;
		  tmpi -= 0.5 * bj;
		  aj = (Re_Data (k1) - Re_Data (k2)) * s60;
		  bj = (Im_Data (k1) - Im_Data (k2)) * s60;
		  Re_Data (k1) = tmpr - bj;
		  Re_Data (k2) = tmpr + bj;
		  Im_Data (k1) = tmpi + aj;
		  Im_Data (k2) = tmpi - aj;
		  kk = k2 + kspan;
	       } while (kk < nn);
	       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 {
		  REALFIX aa, aj, ak, ajm, ajp, akm, akp;
		  REALFIX bb, bj, bk, bjm, bjp, bkm, bkp;
		  int k2, k4;

		  k1 = kk + kspan;
		  k2 = k1 + kspan;
		  k3 = k2 + kspan;
		  k4 = k3 + kspan;
		  akp = Re_Data (k1) + Re_Data (k4);
		  akm = Re_Data (k1) - Re_Data (k4);
		  bkp = Im_Data (k1) + Im_Data (k4);
		  bkm = Im_Data (k1) - Im_Data (k4);
		  ajp = Re_Data (k2) + Re_Data (k3);
		  ajm = Re_Data (k2) - Re_Data (k3);
		  bjp = Im_Data (k2) + Im_Data (k3);
		  bjm = Im_Data (k2) - Im_Data (k3);
		  aa = Re_Data (kk);
		  bb = Im_Data (kk);
		  Re_Data (kk) = aa + akp + ajp;
		  Im_Data (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_Data (k1) = ak - bj;
		  Re_Data (k4) = ak + bj;
		  Im_Data (k1) = bk + aj;
		  Im_Data (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_Data (k2) = ak - bj;
		  Re_Data (k3) = ak + bj;
		  Im_Data (k2) = bk + aj;
		  Im_Data (k3) = bk - aj;
		  kk = k4 + kspan;
	       } while (kk < nn);
	       kk -= nn;
	    } while (kk <= kspan);
	    break;

	  default:
	    k = factor [ii - 1];
	    if (jf != k)
	      {
		 jf = k;
		 s1 = pi2 / (double) jf;
		 c1 = cos (s1);
		 s1 = sin (s1);
		 if (jf > maxFactors)
		   goto Memory_Error;
		 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 {
		  REALFIX aa, ak;
		  REALFIX bb, bk;
		  int k2;

		  aa = ak = Re_Data (kk);
		  bb = bk = Im_Data (kk);

		  k1 = kk;
		  k2 = kk + ispan;
		  j = 1;
		  k1 += kspan;
		  do {
		     k2 -= kspan;
		     Rtmp [j] = Re_Data (k1) + Re_Data (k2);
		     ak += Rtmp [j];
		     Itmp [j] = Im_Data (k1) + Im_Data (k2);
		     bk += Itmp [j];
		     j++;
		     Rtmp [j] = Re_Data (k1) - Re_Data (k2);
		     Itmp [j] = Im_Data (k1) - Im_Data (k2);
		     j++;
		     k1 += kspan;
		  } while (k1 < k2);
		  Re_Data (kk) = ak;
		  Im_Data (kk) = bk;

		  k1 = kk;
		  k2 = kk + ispan;
		  j = 1;
		  do {
		     REALFIX aj = 0.0;
		     REALFIX bj = 0.0;

		     k1 += kspan;
		     k2 -= kspan;
		     jj = j;
		     ak = aa;
		     bk = bb;
		     k = 1;
		     do {
			ak += Rtmp [k] * Cos [jj - 1];
			bk += Itmp [k] * Cos [jj - 1];
			k++;
			aj += Rtmp [k] * Sin [jj - 1];
			bj += Itmp [k] * Sin [jj - 1];
			k++;
			jj += j;
			if (jj > jf)
			  jj -= jf;
		     } while (k < jf);
		     k = jf - j;
		     Re_Data (k1) = ak - bj;
		     Im_Data (k1) = bk + aj;
		     Re_Data (k2) = ak + bj;
		     Im_Data (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 == nFactor)
	   goto Permute_Results;	/* exit infinite loop */
	 kk = jc + 1;
	 do {
	    c2 = 1.0 - cd;
	    s1 = sd;
	    do {
	       c1 = c2;
	       s2 = s1;
	       kk += kspan;
	       do {
		  REALFIX tmp;
		  do {
		     REALFIX ak;
		     ak = Re_Data (kk);
		     Re_Data (kk) = c2 * ak - s2 * Im_Data (kk);
		     Im_Data (kk) = s2 * ak + c2 * Im_Data (kk);
		     kk += ispan;
		  } while (kk <= nt);
		  tmp = s1 * s2;
		  s2 = s1 * c2 + c1 * s2;
		  c2 = c1 * c2 - tmp;
		  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:
   Perm [0] = ns;
   if (kt)
     {
	int k2;

	k = kt + kt + 1;
	if (k > nFactor)
	  k--;
	Perm [k] = jc;
	j = 1;
	do {
	   Perm [j] = Perm [j - 1] / factor [j - 1];
	   Perm [k - 1] = Perm [k] * factor [j - 1];
	   j++;
	   k--;
	} while (j < k);
	k3 = Perm [k];
	kspan = Perm [1];
	kk = jc + 1;
	k2 = kspan + 1;
	j = 1;
	if (nPass != nTotal)
	  {
	     /* permutation for multivariate transform */
	     Permute_Multi:
	     do {
		do {
		   k = kk + jc;
		   do {
		      /* swap
		       * Re_Data (kk) <> Re_Data (k2)
		       * Im_Data (kk) <> Im_Data (k2)
		       */
		      REALFIX tmp;
		      tmp = Re_Data (kk); Re_Data (kk) = Re_Data (k2); Re_Data (k2) = tmp;
		      tmp = Im_Data (kk); Im_Data (kk) = Im_Data (k2); Im_Data (k2) = tmp;
		      kk += inc;
		      k2 += inc;
		   } while (kk < k);
		   kk += (ns - jc);
		   k2 += (ns - jc);
		} while (kk < nt);
		k2 = k2 - nt + kspan;
		kk = kk - nt + jc;
	     } while (k2 < ns);
	     do {
		do {
		   k2 -= Perm [j - 1];
		   j++;
		   k2 = Perm [j] + k2;
		} while (k2 > Perm [j - 1]);
		j = 1;
		do {
		   if (kk < k2)
		     goto Permute_Multi;
		   kk += jc;
		   k2 += kspan;
		} while (k2 < ns);
	     } while (kk < ns);
	  }
	else
	  {
	     /* permutation for single-variate transform (optional code) */
Permute_Single:
	     do {
		/* swap
		 * Re_Data (kk) <> Re_Data (k2)
		 * Im_Data (kk) <> Im_Data (k2)
		 */
		REALFIX t;
		t = Re_Data (kk); Re_Data (kk) = Re_Data (k2); Re_Data (k2) = t;
		t = Im_Data (kk); Im_Data (kk) = Im_Data (k2); Im_Data (k2) = t;
		kk += inc;
		k2 += kspan;
	     } while (k2 < ns);
	     do {
		do {
		   k2 -= Perm [j - 1];
		   j++;
		   k2 = Perm [j] + k2;
		} while (k2 > Perm [j - 1]);
		j = 1;
		do {
		   if (kk < k2)
		     goto Permute_Single;
		   kk += inc;
		   k2 += kspan;
		} while (k2 < ns);
	     } while (kk < ns);
	  }
	jc = k3;
     }

   if ((kt << 1) + 1 >= nFactor)
     FFTRADIX_RETURN(0);
   ispan = Perm [kt];

   /* permutation for square-free factors of n */
   j = nFactor - kt;
   factor [j] = 1;
   do {
      factor [j - 1] *= factor [j];
      j--;
   } while (j != kt);
   nn = factor [kt] - 1;
   kt++;
   if (nn > maxPerm)
     goto Memory_Error;

   j = jj = 0;
   for (;;) {
      int k2;

      k = kt + 1;
      k2 = factor [kt - 1];
      kk = factor [k - 1];
      j++;
      if (j > nn)
	break;				/* exit infinite loop */
      jj += kk;
      while (jj >= k2)
	{
	   jj -= k2;
	   k2 = kk;
	   kk = factor [k++];
	   jj += kk;
	}
      Perm [j - 1] = jj;
   }
   /* determine the permutation cycles of length greater than 1 */
   j = 0;
   for (;;) {
      do {
	 kk = Perm [j++];
      } while (kk < 0);
      if (kk != j)
	{
	   do {
	      k = kk;
	      kk = Perm [k - 1];
	      Perm [k - 1] = -kk;
	   } while (kk != j);
	   k3 = kk;
	}
      else
	{
	   Perm [j - 1] = -j;
	   if (j == nn)
	     break;		/* exit infinite loop */
	}
   }

   maxFactors *= 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 (Perm [j - 1] < 0);
	 jj = jc;
	 do {
	    int k2;

	    if (jj < maxFactors) kspan = jj; else kspan = maxFactors;

	    jj -= kspan;
	    k = Perm [j - 1];
	    kk = jc * k + ii + jj;

	    k1 = kk + kspan;
	    k2 = 0;
	    do {
	       Rtmp [k2] = Re_Data (k1);
	       Itmp [k2] = Im_Data (k1);
	       k2++;
	       k1 -= inc;
	    } while (k1 != kk);

	    do {
	       k1 = kk + kspan;
	       k2 = k1 - jc * (k + Perm [k - 1]);
	       k = -Perm [k - 1];
	       do {
		  Re_Data (k1) = Re_Data (k2);
		  Im_Data (k1) = Im_Data (k2);
		  k1 -= inc;
		  k2 -= inc;
	       } while (k1 != kk);
	       kk = k2;
	    } while (k != j);

	    k1 = kk + kspan;
	    k2 = 0;
	    do {
	       Re_Data (k1) = Rtmp [k2];
	       Im_Data (k1) = Itmp [k2];
	       k2++;
	       k1 -= inc;
	    } while (k1 != kk);
	 } while (jj);
      } while (j != 1);
   }
   FFTRADIX_RETURN(0);			/* exit point here */

   /* alloc or other problem, do some clean-up */
Memory_Error:
   fprintf (stderr, "Error: " FFTRADIXS "() - insufficient memory.\n");
   FFTRADIX_RETURN(-1);
#undef FFTRADIX_RETURN
#undef FFTRADIX_MALLOC
}

/*
 *
 */
int
FFTN (size_t ndim,
      const size_t dims [],
      REAL Re [],
      REAL Im [],
      int iSign,
      REAL scaling)
{
   size_t nTotal;
   int maxFactors, maxPerm;

   /*
    * tally the number of elements in the data array
    * and determine the number of dimensions
    */
   nTotal = 1;
   if (ndim)
     {
	if (dims != NULL)
	  {
	     int i;
	     /* number of dimensions was specified */
	     for (i = 0; i < ndim; i++)
	       {
		  if (dims [i] <= 0) goto Dimension_Error;
		  nTotal *= dims [i];
	       }
	  }
	else
	  nTotal *= ndim;
     }
   else
     {
	int i;
	/* determine # of dimensions from zero-terminated list */
	if (dims == NULL) goto Dimension_Error;
	for (ndim = i = 0; dims [i]; i++)
	  {
	     if (dims [i] <= 0)
	       goto Dimension_Error;
	     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?
    */
   if (dims != NULL)
     {
	int i;
	for (maxFactors = maxPerm = 1, i = 0; i < ndim; i++)
	  {
	     if (dims [i] > maxFactors) maxFactors = dims [i];
	     if (dims [i] > maxPerm) maxPerm = dims [i];
	  }
     }
   else
     {
	maxFactors = maxPerm = nTotal;
     }
#else
   /* use the constants used in the original Fortran code */
   maxFactors = 23;
   maxPerm = 209;
#endif
   /* loop over the dimensions: */
   if (dims != NULL)
     {
	size_t nSpan = 1;
	int i;

	for (i = 0; i < ndim; i++)
	  {
	     int ret;
	     nSpan *= dims [i];
	     ret = FFTRADIX (Re, Im, nTotal, dims [i], nSpan, iSign,
			     maxFactors, maxPerm);
	     /* exit, clean-up already done */
	     if (ret)
	       return ret;
          }
     }
   else
     {
	int ret;
	ret = FFTRADIX (Re, Im, nTotal, nTotal, nTotal, iSign,
			maxFactors, maxPerm);
	/* exit, clean-up already done */
	if (ret)
	  return ret;
     }

   /* Divide through by the normalizing constant: */
   if (scaling && scaling != 1.0)
     {
	int i;

	if (iSign < 0) iSign = -iSign;
	if (scaling < 0.0)
	  scaling = (scaling < -1.0) ? sqrt (nTotal) : nTotal;

	scaling = 1.0 / scaling;	/* multiply is often faster */
	for (i = 0; i < nTotal; i += iSign)
	  {
	     Re_Data (i) *= scaling;
	     Im_Data (i) *= scaling;
	  }
     }
   return 0;

   Dimension_Error:
   fprintf (stderr, "Error: " FFTNS "() - dimension error\n");
   return -1;
}
#endif	/* _FFTN_C */
/*----------------------- end-of-file (C source) -----------------------*/