Why is the transform in Schönhage–Strassen's multiplication algorithm cheap?

The Schönhage–Strassen multiplication algorithm works by turning multiplications of size $N$ into many multiplications of size $lg(N)$ with a number-theoretic transform, and recursing. At least I think that's what it does because there's some other cleverness and I really don't understand it well enough to summarize accurately. It finishes in $O(N \cdot lg(N) \cdot lg(lg(N)))$ time.

A number-theoretic transform is exactly like a Discrete Fourier Transform, except it's done in the finite field $F_{2^N+1}$ of integers modulo $2^N+1$. This makes the operations a lot cheaper, since the Fourier transform has a lot of multiplying by roots of unity, and $F_{2^N+1}$'s roots of unity are all powers of 2 so we can just shift! Also, integers are a lot easier to work with than floating point complex numbers.

Anyways, the thing that confuses me is that $F_{2^N+1}$ is very large. If I give you a random element from $F_{2^N+1}$, it takes $O(N)$ bits to specify it. So adding two elements should take $O(N)$ time. And the DFT does a lot of adding.

Schönhage–Strassen splits the input into $\frac{N}{lg(N)}$ groups with $lg(N)$ bits. These groups are the values of $F_{2^N+1}$ that it will transform. Each pass of the DFT will have $O(\frac{N}{lg(N)})$ additions/subtractions, and there are $O(lg(\frac{N}{lg(N)}))$ passes. So based on addition taking $O(N)$ time it seems like the cost of all those additions should be $O(N \frac{N}{lg(N)} lg(\frac{N}{lg(N)}))$, which is asymptotically the same as $O(N^2)$.

We can do a bit better than that... because the values start out so small, the additions are quite sparse. The first pass' additions really only cost $O(lg(N))$ each, and the second pass' cost $2^1 O(lg(N))$ each, and the i'th pass' cost $O(min(N, 2^i \cdot lg(N)))$ each, but that still all totals out to a terrible $\frac{N^2}{lg(N)}$.

How is Schönhage–Strassen making the additions cheap? How much do they cost overall?

Is it because the algorithm actually uses $F_{N+1}$ (with $N$ guaranteed to be a power of 2)? There's enough stacked $2^{2^{n}}$ and german in the paper that I'm really not sure. On the other hand, I don't think that guarantees enough roots of unity for things to work out.

• Have you read the Wikipedia article? At the very least, it's not in German. According to Wikipedia, the size is reduced from $N$ to $n \approx \sqrt{N}$ at each step, which works out since then $n^2 = \Theta(N)$. – Yuval Filmus Jul 11 '14 at 3:00
• @Yuval Where do you see $n \approx \sqrt{N}$ in the article? There's a lot of "primitive root", and an extra optimization trick involving square roots, but otherwise I don't see it. I do see "Split each input number into vectors X and Y of $2^k$ parts each, where $2^k$ divides N.", but no definition of $k$. – Craig Gidney Jul 11 '14 at 3:09
• @Yuval That does seem to help... but there's still $O(\sqrt{N})$ additions per pass with each costing $O(N)$. So a pass costs $O(N^{\frac{3}{2}})$, which is better but not good enough. – Craig Gidney Jul 11 '14 at 3:21
• The setting of $k$ appears in the sentence "The optimal number of pieces to divide the input into is proportional to $\sqrt{N}$." The value of $n$ is then $2^{\Theta(\sqrt{N})}+1$, so each addition costs $O(\sqrt{N})$. – Yuval Filmus Jul 11 '14 at 4:00

According to the Wikipedia article, at each step the length of the integers is reduced from $N$ to (roughly) $\sqrt{N}$, and there are (roughly) $\sqrt{N}$ of them, and so the additions only cost $O(N)$. There is a detailed analysis of the running time in the final paragraph of the linked section, copied here in case it changes:

In the recursive step, the observation is used that:

• Each element of the input vectors has at most $N/2^k$ bits;
• The product of any two input vector elements has at most $2N/2^k$ bits;
• Each element of the convolution is the sum of at most $2^k$ such products, and so cannot exceed $2N/2^k + k$ bits.

Here $N$ is current input length and $k = \Theta(\sqrt{N})$. Arithmetic is done modulo $2^n+1$, where $n$ is some multiple of $2^k$ which is larger than $2N/2^k + k$; note that $n = \Theta(\sqrt{N})$.

• I don't understand why each element of the convolution can't exceed $\frac{2N}{2^k} + k$ bits. I understand that if we just added them together without any shifting that they would be no larger than that, but aren't we doing arithmetic in $F_{2^N+1}$ and scaling things by huge powers of 2? The $\sqrt{N}$'th element of the convolution is $\sum_{i=0}^{\sqrt{N}} x_i 2^{i \sqrt{N}}$, stepping up the power by $\sqrt{N}$ for each element, so won't it fail to benefit from overlapping their bits at all and require $N$ bits? – Craig Gidney Jul 11 '14 at 4:22
• Are the additions being doing in $F_{2^N+1}$? The article says "it's necessary to use a smaller N for recursive multiplications", but it's not clear if that step down is done for both the transform and the pointwise multiplications or only for the pointwise multiplications. – Craig Gidney Jul 11 '14 at 4:44
• The Wikipedia article mixes up $N$ and $n$; hopefully that's fixed in my revised answer. Arithmetic is done modulo $2^n+1$, where $n = \Theta(\sqrt{N})$. The resulting vector can then be combined to give a result modulo $2^N+1$. – Yuval Filmus Jul 11 '14 at 4:52
• Reiterating, the trick is to perform arithmetic modulo a smaller modulus, whose size is compatible with the size of the pieces. The $N$-bit integers are cut into roughly $\sqrt{N}$ pieces of length $\sqrt{N}$ bits. – Yuval Filmus Jul 11 '14 at 4:54
• Ah, that explains it then. That would mean the $lg(n)$ version is also quite fast, since it has $lg(n)$ sized operations. If you could explicitly say that in your answer, that'd be great. – Craig Gidney Jul 11 '14 at 4:57

There's an excellent explanation of exactly what's going on, including how the size of the field goes down as you progress, in GMP's documentation: 15.1.6 FFT Multiplication.

A more theoretically in-depth overview starts on page 55 of Modern Computer Arithmetic (pdf) (well... the browser says page 71 but the page itself says 55). It explains the algorithm's correctness and breaks down the complexity analysis.