Schönhage–Strassen multiplication

Schönhage and Strassen [SS1971] multiply integers with a discrete Fourier transform over the ring \mathbb{Z}/(2^N + 1)\mathbb{Z}. In this ring 2 has order 2N, so the roots of unity needed by the FFT can be chosen as powers of two, and multiplying by them is a shift with wraparound, costing O(N) bit operations instead of a multiplication. The pointwise products are products in the same kind of ring and can be computed with the same method, which gives O(n \log n \log \log n) bit operations for n-bit integers. Fürer [Fur2009] and Harvey and van der Hoeven [HvdH2021] improved the asymptotic bound to O(n \log n), but the Schönhage–Strassen algorithm remains the standard exact integer FFT, as in GMP [GKZ2007] and in FLINT's fft module, written by William Hart.

The transforms of fft.h are truncated (Truncated FFTs) and, for large lengths, computed with the matrix Fourier algorithm (The matrix Fourier algorithm).

The ring#

Let n = 2^{\mathit{depth}}, let w \ge 1, and let N = nw and p = 2^N + 1. Put \theta = 2^w \in \mathbb{Z}/p\mathbb{Z}.

The lemma says that \theta is a principal root of unity, which is what the inversion formula needs in a ring that is not a field [Pol1971]: for \hat a_k = \sum_j a_j \theta^{jk},

and 2n is invertible. Hence pointwise products of transforms are the transforms of cyclic convolutions of length 2n, and the inverse transform recovers these after a division by 2n, which is a shift.

Multiplication by \theta^i = 2^{iw} is a shift by iw bits. Since 2^{N + s} \equiv -2^s, the bits shifted out at the top re-enter at the bottom with a negative sign. FLINT splits the shift into whole limbs, realised by reading or writing at an offset with the wrapped part subtracted, and a remaining shift by fewer than \mathtt{FLINT\_BITS} bits. Residues are kept in N/\mathtt{FLINT\_BITS} + 1 limbs in two's complement, and sums and differences are not reduced: the top limb absorbs the carries, and mpn_normmod_2expp1 reduces only where a reduced value is needed. A butterfly (a, b) \mapsto (a + b, 2^{iw}(a - b)) is thus a combined addition and subtraction (flint_mpn_sumdiff_n) followed by a shift: a few passes over N bits. With these butterflies the transforms are radix-2 Cooley–Tukey FFTs (Radix-2 Cooley–Tukey).

The \sqrt 2 trick#

A transform of length 2n with root 2^w needs a ring of N = nw \ge n bits. The following identity, also used by GMP [GKZ2007], gives a root of unity of order 4n in the same ring.

Even powers \varphi^{2j} = 2^{jw} are shifts. For odd i and odd w, iw is odd and \varphi^i = \sigma^{iw} = \sigma \cdot 2^{(iw - 1)/2}. Since \sigma = 2^{N/4} (2^{N/2} - 1), multiplying by \varphi^i costs a shift by (iw - 1)/2 + N/4 bits, a shift by N/2 bits and a subtraction; this is fft_adjust_sqrt2 and fft_butterfly_sqrt2. For even w, \varphi = 2^{w/2} and nothing special is needed. In a transform of length 4n by decimation in frequency, only the first layer uses odd powers of \varphi; the two halves are transforms of length 2n with root \theta. The trick therefore doubles the transform length for a given ring, or halves the ring for a given number of points, at the cost of a slightly more expensive first layer. All _sqrt2 functions of fft.h use it, with N a multiple of \mathtt{FLINT\_BITS}, hence of 4.

Integer multiplication#

To multiply x and y, FLINT cuts them into chunks of b bits, x = \sum_{i < j_1} x_i 2^{bi} and y = \sum_{i < j_2} y_i 2^{bi} with 0 \le x_i, y_i < 2^b (fft_split_bits). The product xy is the value at 2^b of the polynomial \sum_k c_k t^k with c_k = \sum_{i + j = k} x_i y_j, which is computed by a convolution of length 4n modulo 2^N + 1 using the \sqrt 2 trick.

The bound 2^{\mathit{depth} + 1 + 2b} \le 2^N means that almost half of each transformed coefficient carries data: N \approx 2b. The product has about 4n \cdot b \approx 2nN bits, so each transform handles about twice the size of the product. mul_mfa_truncate_sqrt2 performs the same steps with the matrix Fourier algorithm, and fft_convolution performs steps 3 to 5 on given coefficient vectors, also with signed coefficients, which are then read from their residues in the symmetric range (-p/2, p/2). When only a window of the coefficients is wanted, as in a middle product, fft_convolution_choose_mid shortens the cyclic convolution so that the wrap-around misses the window (Partial products), as the small-prime FFT does (Windows and wrap-around).

Multiplication modulo 2^N + 1#

The pointwise products in step 4 are products modulo 2^N + 1. Up to FFT_MULMOD_2EXPP1_CUTOFF limbs (128 limbs, 8192 bits), flint_mpn_mulmod_2expp1_basecase computes the full product with flint_mpn_mul_n and subtracts its high half from its low half. Beyond that, fft_mulmod_2expp1 applies the method of Schönhage and Strassen once more, now with a negacyclic convolution, which reduces modulo 2^N + 1 for free.

Write N = 2n' b' with n' = 2^{d'} and cut x and y, both less than 2^N, into 2n' pieces x_i, y_i of b' bits. Since (2^{b'})^{2n'} = 2^N \equiv -1, the product xy modulo 2^N + 1 is the value at 2^{b'} of X(t) Y(t) modulo t^{2n'} + 1, whose coefficients are

so that |c_k| < 2n' \, 2^{2b'}. The negacyclic convolution is a cyclic convolution of the weighted vectors (\varphi'^i x_i) and (\varphi'^i y_i), where \varphi' is a root of unity of order 4n', followed by multiplication of coefficient i by \varphi'^{-i} [Nus1980; CF1994]. FLINT computes it modulo p' = 2^{N'} + 1 with N' = 2b' = n' w' and \varphi' = \sqrt 2^{\,w'} (fft_negacyclic, ifft_negacyclic). This ring is too small for the coefficients by d' + 2 bits; the missing information comes from a second, cheap convolution.

The coefficients satisfy |c_k| < 2n' \, 2^{N'} < (B/2 - 1) p' as long as 4n' < B - 2, which holds for any feasible transform. The residues modulo B are the negacyclic convolution of the lowest limbs of the pieces in wrapping word arithmetic, fft_naive_convolution_1, which costs 4n'^2 word multiplications, negligible next to the rest. Without this trick the inner ring would need 2b' + d' + 2 bits, which is not of the form n' w' with a multiple of \mathtt{FLINT\_BITS}. Finally, xy \bmod (2^N + 1) is assembled from the c_k as in step 6 of Algorithm 1, with a correction for the negative ones.

fft_adjust_limbs makes step 3 well defined: it rounds the limb count up so that N is a multiple of n'^2, making w' an integer, and b' a multiple of \mathtt{FLINT\_BITS}. The pointwise products of step 4 are not split further, so fft_mulmod_2expp1 has one level of recursion; its inner products of 2b' bits are full products by flint_mpn_mul_n. In builds with fft_small, the convolutions of fmpz_poly and fmpzi use instead sd_fft_mpn_mulmod_2expp1, a negacyclic convolution modulo word-size primes in floating point, whenever fft_adjust_limbs_sd_fft finds a suitable ring; it is about twice as fast for N \ge 16384. See The negacyclic engine of fft_small and, for the algorithm, Algorithm.

Choice of parameters#

flint_mpn_mul_fft_main chooses (\mathit{depth}, w) in two steps.

1. It tries (\mathit{depth}, w) in the order (6, 1), (6, 2), (7, 1), (7, 2), \ldots and stops at the first pair for which chunks of b = \lfloor (nw - \mathit{depth} - 1)/2 \rfloor bits fit, j_1 + j_2 - 1 \le 4n.
2. If \mathit{depth} < 11, it shortens the transform by the factor 2^k given by the table entry FFT_TAB[depth - 6][w - 1] and multiplies w by 4^k, which keeps the capacity 4n \cdot b roughly the same with 2^k times fewer, larger coefficients. It then lowers w in steps of \max(1, 2^{6 - \mathit{depth}}), keeping nw a multiple of 64, as long as the product still fits, and calls mul_truncate_sqrt2. Otherwise, if j_1 + j_2 - 1 \le 3n, it halves the transform length and triples w, and then calls mul_mfa_truncate_sqrt2.

For small products, few large coefficients win because the pointwise products, by Toom–Cook and basecase multiplication, are faster per bit than additional transform layers; for large products, the transform must stay long so that the coefficients and the pointwise products remain small. The tables are tuned per processor with src/fft/tune/tune-fft.c. For \mathit{depth} = 6, \ldots, 10, each with w = 1, 2:

With the generic table, balanced products of \ell limbs each get the following parameters ("first fit" is the pair found in step 1, "truncate" and "MFA" stand for mul_truncate_sqrt2 and mul_mfa_truncate_sqrt2):

The pointwise products of fft_mulmod_2expp1 use d' = \lfloor \lceil \log_2 N \rceil / 2 \rfloor - k with k from MULMOD_TAB, indexed by \lceil \log_2 N \rceil - 12 and clamped at the last entry. The values lie between 1 and 4; the generic table has 19 entries, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 1, 1. The pieces then have between about \sqrt N and 8 \sqrt N bits, as in the last column above. In builds without fft_small, flint_mpn_mul switches from GMP to flint_mpn_mul_fft_main at FLINT_FFT_MUL_THRESHOLD = 32000 limbs.

Complexity#

In this section n is the bit size of the operands. A transform of length L with N-bit coefficients costs O(L \log L) butterflies of O(N) bit operations. With LN \approx 4n, the transforms cost O(n \log n) bit operations, and the L pointwise products cost L \cdot M(N), where M(N) is the cost of a multiplication modulo 2^N + 1. Choosing N \approx \sqrt n and applying the method recursively to the pointwise products gives the bound O(n \log n \log \log n) of [SS1971], since each of the O(\log \log n) levels costs O(n \log n). FLINT stops after one negacyclic level and uses the best available multiplication for the remaining products of a few thousand bits. In practice the cost per bit grows slowly with n and is dominated by the memory traffic of the transforms, which the matrix Fourier algorithm keeps in cache. FFTs modulo word-size primes in floating point are faster on current processors (FFTs modulo word-size primes in floating point); FFTs in FLINT compares the methods.