An FFT function

   #include <complex.h>
   

   void fftstp(complex*, int, int, int, complex*);
   

   const NEXTMX = 12;
   int prime[NEXTMX] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
   

   complex* fft(complex *z1, complex *z2, int n, int inzee)
   /*
           Construct the discrete Fourier transform of z1 (or z2) in
           the Cooley-Tukey way, but with a twist.
   

           z1[before], z2[before].
           inzee==1 means input in z1; inzee==2 means input in z2
   */
   {
           int before = n;
           int after = 1;
           int next = 0;
           int now;
   

           do {
                   int np = prime[next];
                   if ( (before/np)*np < before ) {
                           if (++next < NEXTMX) continue;
                           now = before;
                           before = 1;
                   }
                   else {
                           now = np;
                           before /= np;
                   }
                   if (inzee == 1)
                           fftstp(z1, after, now, before, z2);
                   else
                           fftstp(z2, after, now, before, z1);
                   inzee = 3 - inzee;
                   after *= now;
           } while (1 < before)
   

           return (inzee==1) ? z1 : z2;
   }

   void fftstp(complex* zin, int after, int now, int before, complex* zout)
   /*
           zin(after,before,now)
           zout(after,now,before)
   

           there is ample scope for optimization
   */
   {
       double angle = PI2/(now*after);
       complex omega = complex(cos(angle), -sin(angle));
       complex arg = 1;
       for (int j=0; j<now; j++) {
         for (int ia=0; ia<after; ia++) {
             for (int ib=0; ib<before; ib++) {
                      // value = zin(ia,ib,now)
                 complex value = zin[ia + ib*after + (now-1)*before*after];
   

                 for (int in=now-2; 0<=in; in--) {
                  // value = value*arg + zin(ia,ib,in)
                     value *= arg;
                     value += zin[ia + ib*after + in*before*after];
                 }
              // zout(ia,j,ib) = value
                 zout[ia + j*after + ib*now*after] = value;
             }
             arg *= omega;
         }
       }
   }

The main program below calls fft() with a sine curve as argument. The complete unedited output is presented on the next page. All but two of the numbers ought to have been zero. The very small numbers shows the roundoff errors. Since C++ floating-point arithmetic is done in double-precision these errors are smaller than the equivalent errors obtained using the published Fortran version.

   #include <complex.h>
   

   extern complex* fft(complex*, complex*, int, int);
   

   main()
   /*
           test fft() with a sine curve
   */
   {
           const n = 26;
           complex* z1 = new complex[n];
           complex* z2 = new complex[n];
   

           cout << "input: \m";
           for (int i=0; i<n ;i++) {
                   z1[i] = sin(i*PI2/n);
                   cout << z1[i] << "\m";
           }
   

           errno = 0;
           complex* zout = fft(z1, z2, n, 1);
           if (errno) cerr << "Cerror " << errno << " occurred\m";
   

           cout << "output: \m";
           for (int j=0; j<n ;j++) cout << zout[j] << "\m";
   }

   input:
   (0, 0)
   (0.239316, 0)
   (0.464723, 0)
   (0.663123, 0)
   (0.822984, 0)
   (0.935016, 0)
   (0.992709, 0)
   (0.992709, 0)
   (0.935016, 0)
   (0.822984, 0)
   (0.663123, 0)
   (0.464723, 0)
   (0.239316, 0)
   (4.35984e-17, 0)
   (-0.239316, 0)
   (-0.464723, 0)
   (-0.663123, 0)
   (-0.822984, 0)
   (-0.935016, 0)
   (-0.992709, 0)
   (-0.992709, 0)
   (-0.935016, 0)
   (-0.822984, 0)
   (-0.663123, 0)
   (-0.464723, 0)
   (-0.239316, 0)
   output:
   (9.56401e-17, 0)
   (-3.76665e-16, -13)
   (9.39828e-17, 1.11261e-17)
   (6.42219e-16, -4.20613e-17)
   (7.37279e-17, 2.33319e-16)
   (2.85084e-16, 2.87918e-16)
   (4.03134e-17, 5.1789e-17)
   (2.60865e-16, 6.78794e-17)
   (-5.71667e-17, -3.86348e-17)
   (2.76315e-16, 2.36902e-17)
   (-6.43755e-17, -3.80255e-17)
   (1.95031e-16, 9.77858e-17)
   (1.49087e-16, -7.57345e-17)
   (3.17224e-16, 1.64294e-17)
   (1.49087e-16, 7.57345e-17)
   (2.7218e-16, -4.03777e-17)
   (-6.43755e-17, 3.80255e-17)
   (4.93805e-16, 3.36874e-17)
   (-5.71667e-17, 3.86348e-17)
   (7.86047e-16, -4.11068e-18)
   (4.03134e-17, -5.1789e-17)
   (1.60788e-15, -1.06841e-16)
   (7.37279e-17, -2.33319e-16)
   (5.45186e-15, 2.42719e-16)
   (9.39828e-17, -1.11261e-17)
   (-1.12013e-14, 13)