FFT via ring morphisms.

gistfile1.txt

Pick a 2n-element periodic sequence of complex numbers. This is isomorphic to a polynomial p in
  C[x] / (x^2n - 1)
(x^2n = 1 = x^0, which encodes the 2n-periodicity).

Computing the 2n-element DFT of the original sequence is just evaluating p at the roots of unity
(\omega_{2n}^0, \omega_{2n}^1, ..., \omega_{2n}^{2n+1}).

Now, (x^2n - 1) = (x^n - 1) (x^n + 1). This implies that the given p mod (x^2n - 1) is uniquely
determined by the remainders p mod (x^n - 1) and p mod (x^n + 1) [via CRT]. That is, there is
an isomorphism between C[x] / (x^2n - 1) and (C[x] / (x^n - 1)) x (C[x] / (x^n + 1)).

How do we compute these remainders? By substitution! Given some p mod (x^2n - 1), we have

  p = p_0 + p_1 x + p_2 x^2 + ... + p_{2n-1} x^{2n-1}
  
Computing p mod (x^n - 1) just means saying that (x^n - 1) = 0, i.e. x^n = 1. That is, in the
second half of the above sum, we can factor out the x^n and replace it by 1; collect terms and
we get

  p mod (x^n - 1) = (p_0 + p_n) + (p_1 + p_{n+1}) x + ... + (p_{n-1} + p_{2n-1}) x^{n-1}
  
Likewise, p mod (x^n + 1) is very similar, except now we have x^n = -1, so we get a sign flip:

  p mod (x^n + 1) = (p_0 - p_n) + (p_1 - p_{n+1}) x + ... + (p_{n-1} - p_{2n-1}) x^{n-1}

And if we wanted to, then for any power-of-2 initial n, we could factor any

  (x^2n - r^2) = (x^n - r) (x^n + r)
  
and then determine

  p mod (x^n - r) = (p_0 + r * p_n) + (p_1 + r * p_{n+1}) + ... + (p_{n-1} + r * p_{2n-1}) x^{n-1}
  p mod (x^n + r) = (p_0 - r * p_n) + (p_1 - r * p_{n+1}) + ... + (p_{n-1} - r * p_{2n-1}) x^{n-1}
  
We already have a fully functional power-of-2 size FFT algorithm (Gauss's) at this point. Just
keep splitting the polynomial by ever-smaller roots of unity, until we're down to linear factors.
Note that a linear factor

  p mod (x - r) = k
  
just is a fancy way of saying that p(r) = k (polynomial remainder theorem); once we've split
everything down to linear factors at the roots of unity, we've calculated p's value at these
points, which was our original task. All but the first few steps end up with non-trivial r's
though. This is the equivalent of the twiddle factors, although the twiddles in more typical
FFT algorithms do something different, see below.

However, while this is _a_ FFT algorithm, it's not the ones you usually see, and it's desirable
to get to the usual recursive decomposition where you compute sub-FFTs then merge them. How do
we get there?
  
If we look at our original split, the "p mod (x^n - 1)" part is all good; in the same shape we
started with, ready to go for another round. The problem is "p mod (x^n + 1)", which lives in
a different ring. Can we fix that?

Why, yes! We scale the input argumnet to get a slightly different polynomial in our original ring.
Namely, we go from p(x) to p'(x) = p(\omega_{2n} x). Now with our original

  p = p_0 + p_1 x + ... + p_{n-1} x^{n-1}    [mod (x^n + 1)]

we have

  p' = p_0 + p_1 \omega_{2n} x + ... + p_{n-1} \omega_{2n}^{n-1} x^{n-1}    [mod (x^n - 1)]
  
where the powers of \omega_{2n} are just from the scaling of x. The modulus changes because if
x^n = -1 then (\omega_{2n} x)^n = (\omega_{2n}^n) x^n = -1 * -1 = 1.

...and that's our twiddle factors right there. The even part stays as it is. The odd part
starts out in a different ring; but if we multiply the coefficients by the twiddles, we're
back to a ring with the same structure as the even part (so we can call the FFT recursively).

The recursive decomposition thus obtained is decimation-in-frequency, which is more natural in
this setting.

Why bother doing it like this? One of the nice things about this approach is that in this form,
some algorithmic variants are a lot cleaner to understand. What I've shown so far is:

1. Gauss-style FFT: deal with arbitrary roots r in the split.
2. "Regular" FFT style: factor (x^2n - 1) -> (x^n - 1) x (x^n + 1), then map the second part
   to (x^n - 1) by the x |-> \omega_{2n} x change of coordinates, recurse on both halves.
   
variant 2 has all the multiplies in the change of coordinates; the two remainder computations
are just additions and subtractions.

However, the next step in the decomposition of a Gauss-style FFT has r=+-i, which is still
cheap, since the multiplications by +-i are "free" in terms of arithmetic performed on real
numbers (the "p_k + r * p_{n+k}" calculations still boil down to just additions and subtractions).
Therefore, we can stay in that form for one more stage and only do the change of coordinates
after; this yields the split-radix FFT:

Factor (x^4n - 1) -> (x^2n - 1) * (x^2n + 1) -> (x^2n - 1) * (x - i) * (x + i) and then do
the change of coordinates on the two latter parts, at which point we're back to a ring
factorization C[x]/(x^2n - 1)  x  C[x]/(x^n - 1)  x  C[x]/(x^n - 1)  and can recurse. The win
here is that we save a round of twiddling in the middle step.

----

I got this view of the FFT from the first few sections of Dan Bernstein's "The tangent
FFT" [1], which elaborates on his much terser presentation in section 7 of "Multidigit
multiplication for Mathematicians" [2].

It's not a common presentation (hence my going into details here) but I agree with
Bernstein that it is extremely clean and powerful.

[1] http://cr.yp.to/arith/tangentfft-20070919.pdf
[2] http://cr.yp.to/papers/m3.pdf