Listing 1 Tuning an FFT by simplifying the use of trig functions

#include <math.h>
#define SWAP(a,b) tempr=a;a=b;b=tempr
void four1(float *data, int *nn, int *isign)
{                               /* altered for consistency
                           * with original FORTRAN */
   /* Press, Flannery, Teukolsky, Vettering "Numerical
    * Recipes in C" tuned up ; Code works only when *nn is
    * a power of 2  */
   int n, mmax, m, j, i;
   double wtemp, wr, wpr, wpi, wi, theta, wpin;
   float tempr, tempi, datar, datai;
   n = *nn * 2;
   j = 0;
   for (i = 0; i < n; i += 2) {
      if (j > i) {             /* could use j>i+1 to help
                            * compiler analysis */
         SWAP(data[j], data[i]);
         SWAP(data[j + 1], data[i + 1]);
      }
      m = *nn;
      while (m >= 2 && j >= m) {
         j -= m;
         m >>= 1;
      }
      j+=m;
   }
   theta = 3.141592653589795 * .5;
   if (*isign < 0)
      theta = -theta;
   wpin = 0;                   /* sin(+-PI) */
   for (mmax = 2; n > mmax; mmax *= 2) {
      wpi = wpin;
      wpin = sin(theta);
      wpr = 1 - wpin * wpin - wpin * wpin;    /* cos(theta*2) */
      theta *= .5;
      wr = 1;
      wi = 0;
      for (m = 0; m < mmax; m += 2) {
         for (i = m; i < n; i += mmax * 2) {
            j = i + mmax;
            /* mixed precision not significantly more
             * accurate here; if removing float casts,
             * tempr and tempi should be double */
            tempr = (float) wr*data[j] - (float) wi *data[j + 1];
            tempi = (float) wr *data[j + 1] + (float) wi *data[j];
            /* don't expect compiler to analyze j >
             * i+1; original code stored data[j..]
             * first and avoided using data[ri] */
            data[i] = (datar= data[i]) + tempr;
            data[i + 1] = (datai = data[i + 1]) + tempi;
            data[j] = datar - tempr;
            data[j + 1] = datai - tempi;
         }
         wr= (wtemp = wr) * wpr - wi * wpi;
         wi = wtemp * wpi + wi * wpr;
      }
   }
}
/* End of File */