Complexity of computing the discrete Fourier transform?

Question

What is the complexity (on the standard integer RAM) of computing the standard discrete Fourier transform of a vector of $n$ integers?

The classical algorithm for fast Fourier transforms, inappropriately^[1] attributed to Cooley and Tukey, is usually described as running in $O(n \log n)$ time. But most of the arithmetic operations executed in this algorithm start with complex $n$th roots of unity, which are (for most $n$) irrational, so exact evaluation in constant time is not reasonable. The same issue arises with the naive $O(n^2)$-time algorithm (multiplying by a Vandermonde matrix of complex roots of unity).

It's not even clear how to represent the output of the DFT exactly (in any useful form). In other words, it's not clear that computing DFTs is actually possible!

So suppose we only need $b$ bits of precision in each output value. What's the complexity of computing the discrete Fourier transform, as a function of $n$ and $b$? (For concreteness, feel free to assume $n$ is a power of $2$.)

Or does every instance of "FFT" in the literature actually mean "fast number-theoretic transform"?^[2]

See my related questions on the complexity of Gaussian elimination and Euclidean shortest paths.

_{^[1] It should really be called (some prefix of) the Gauss-Runge-König-Yates-Stumpf-Danielson-Lánczos-Cooley-Tukey algorithm.}

_{^[2] And if so, why do most textbooks describe only the complex-number algorithm?}

Answer 1

This answer is a variant of the analysis of the first algorithm ("Methode A") by Schönhage and Strassen for multiplication of long integers.

Assume we want to compute an FFT of length $K = 2^k$. Scale your input such that all values are smaller than 1. Let us first assume that we compute with $m$-bit fixed point arithmetic ($m$ bits after the binary point). Let $\delta = 2^{1/2 -m}$ be the ("complex") unit of least position. Let $\omega = \exp(2\pi i/K)$.

1) One can compute approximations $\omega_j'$ such that $|\omega_j' - \omega^j| \le (2k-1)\delta$ for all $0 \le j \le K-1$. This can be done in time $O(K M(m))$ where $M(m)$ is the time needed to multiply $m$-bit numbers. (see Knuth Vol. 2, 3rd ed., page 309).

If standard integer RAM means logarithmic cost, then $M(m) = O(m \log m)$. If standard integer RAM means word RAM, then $M(m) = O(m)$. (Schönhage and Strassen show in "Methode A" how to reduce in linear time the multiplication of $m$-bit numbers to $m$ multiplication of $O(\log m)$ bit numbers. The latter can be done at unit costs.)

2) The classical Cooley-Tukey FFT computes operations of the form $a = b + \omega^j c$. We use $m$-bit fixed point arithmetic, these opertions become $a' = truncate(b' + \omega_j' c')$. If we know $b'$ and $c'$ up to an error of $\epsilon$, we get $a'$ up to an error of $2\epsilon + 2k\delta$.

3) Using induction, it is easy to see that we get the final result with error $(2^k - 1) \cdot 2k\delta$. To get precision $b$ in the end, $m \ge k + \log k + b + O(1)$.

4) Thus the final running time is $O(K k M(k+b))$.

This should also work with floating point numbers: 1) can still be done with fixed point arithmetic, 2) is also true for floating point numbers.

---

In fixed point arithmetic, I think, it can even be done faster. First we reduce the computation of the FFT to the multiplication of polynomials using Bluestein's trick. The length of the coefficients needed to get the desired precision should be $O(k + b)$. Then we reduce the multiplication of polynomials to the multiplication of long integers. (Append the coefficients to a long number and separate them by blocks of zero of length $O(k+b)$.) The length of the integers is $O(K(k+b))$.

Answer 2

This is not a complete answer, but I can point you to some relevant papers and also partially explain why it's not so easy to extract an answer to your specific question from the literature.

Let me start by asking, why do you want to know the answer to this question? Typically, the people who find themselves caring about this sort of issue are those faced with actually implementing a high-performance FFT for a practical application. Such people care less about asymptotic complexity in some idealized computational model than about maximizing performance under their particular hardware and software constraints. For example, the developers of the Fastest Fourier Transform in the West write in their paper:

The best choice depends upon hardware details like the number of registers, latency and throughput of instructions, size and associativity of caches, structure of the processor pipeline, etc.

These are issues that theorists typically don't want to sully their hands with, but they are of great importance in actual implementations. If a theorist declares, "I've figured out the absolute best asymptotic bit complexity in the RAM model," the practitioner might say, "That's nice," but may find such a theoretical result useless for his or her purposes.

Having said that, I think that your best bet is to look at the numerical analysis literature. For example, Tasche and Zeuner have taken a close look at the numerical stability of the FFT algorithm. This may still not be exactly what you want, because the general consensus among practitioners seems to be that to achieve a given amount of numerical precision, the best practical approach is to precompute certain numbers called "twiddle factors" to high accuracy. If you're doing only one FFT, then this is not going to be the fastest approach because you don't get to amortize the cost of your one-time precomputation over a large number of FFT computations. Still, their analysis of the worst-case roundoff error should still be relevant to your question.