A recent slogan is “A program must not only work but must appear to work”. The above program works but alas it does not appear to me to work. There are two difficulties in understanding the short program:

- The mysterious code: “{int m=n/2; while(m<=j) {j -= m; m /= 2;} j += m;}” does not fall into any familiar pattern.
- Even a complete understanding of what this code does, sheds little light on its purpose.
- The original question remains of understanding why the floating point calculations performed produce the required transform. Here is an explanation that is not easy to follow.

The routine “fft” takes a reference to an array of complex numbers and returns when it has replaced those numbers with their discrete Fourier transform. The size must be a power of two. This proof is for the particular size, 1024, of the array. (If your array size is different you can modify and run my code.)

One prerequisite is a little theory on vector spaces, in particular a 1024 dimensional vector space over the field of complex numbers. The 1024 complex numbers that are the inputs to fft are the components of a vector in the space and so are the returned values. The transformation between input in output is linear.

I provide here a top level view of the proof here with links to elaborations of the steps. This is the only glue that holds the proof together.

- I divide the identifiers for numbers into three classes
- Control
- The fixed identifiers (declared “int”)
- bControl
- The complex identifiers “b” & “w”.
- Linear
- These are the complex numbers, (type “c”), aside from b & w.

**values**as the routine runs. - Each Control and bControl identifier takes on a sequence of values during the execution of the routine and these values do not depend on the complex inputs to the routine.
- There are no references to uninitialized variables.
- Array subscript values are always in bounds.
- In summary so far: The order of execution of the statements and expressions of the routine is deterministic. Each control construct, “if”, “for” and “while” is driven by Control values and the sequence of control values does not depend on the complex inputs (a[]).
- The Linear values, several per complex identifier, are each a linear combination of the routine inputs.
- The output values are each a linear combination of the inputs.
- Two linear transformations that agree on each member of some set of basis elements, agree everywhere.
- The routine produces the correct values for each call when the input array is all 0 except for one element with value 1.
- The above set of input possibilities constitute a basis for the space of inputs for the routine.
- The routine is always correct.

This is a multiprocessor version for a 2D transform. See these amusing and useful notes on 2D Fourier transforms. Transformed images