Schönhage–Strassen Algorithm

The Schönhage–Strassen algorithm is an asymptotically fast multiplication algorithm for large integers. It was developed by Arnold Schönhage and Volker Strassen in 1971.A. Schönhage and V. Strassen,
Schnelle Multiplikation großer Zahlen
, ''Computing'' 7 (1971), pp. 281–292. The run-time bit complexity is, in big O notation, $O(n \cdot \log n \cdot \log \log n)$ for two ''n''-digit numbers. The algorithm uses recursive fast Fourier transforms in rings with 2^''n''+1 elements, a specific type of number theoretic transform.

The Schönhage–Strassen algorithm was the asymptotically fastest multiplication method known from 1971 until 2007, when a new method, Fürer's algorithm, was announced with lower asymptotic complexity, and is used in the Basic Polynomial Algebra Subprograms (BPAS) open source library. The algorithm does not adapt for polynomials over finite fields though. The current best multiplication algorithm in terms of asymptotic complexity is by David Harvey and Joris van der Hoeven, which gives $O(n\log n)$ complexity. However, this improvement is too insignificant to be seen in integer sizes seen in practice (see Galactic algorithms) and implementations do not exist.

In practice the Schönhage–Strassen algorithm starts to outperform older methods such as Karatsuba and Toom–Cook multiplication for numbers beyond 2²¹⁵ to 2²¹⁷ (10,000 to 40,000 decimal digits).Overview of Magma V2.9 Features, arithmetic section
: Discusses practical crossover points between various algorithms. The GNU Multi-Precision Library uses it for values of at least 1728 to 7808 64-bit words (33,000 to 150,000 decimal digits), depending on architecture. There is a Java implementation of Schönhage–Strassen which uses it above 74,000 decimal digits.

Applications of the Schönhage–Strassen algorithm include mathematical empiricism, such as the Great Internet Mersenne Prime Search and computing approximations of ''π'', as well as practical applications such as Kronecker substitution, in which multiplication of polynomials with integer coefficients can be efficiently reduced to large integer multiplication; this is used in practice by GMP-ECM for Lenstra elliptic curve factorization.

Implementation details

This section is based primarily on an overview of the method by Crandall and Pomerance in their ''Prime Numbers: A Computational Perspective''.R. Crandall & C. Pomerance. ''Prime Numbers – A Computational Perspective''. Second Edition, Springer, 2005. Section 9.5.6: Schönhage method, p. 502. This variant differs somewhat from Schönhage's original method in that it exploits the discrete weighted transform to perform negacyclic convolutions more efficiently. Another source for detailed information is Knuth's ''The Art of Computer Programming''.Donald E. Knuth, The Art of Computer Programming, Volume 2: Seminumerical Algorithms (3rd Edition), 1997. Addison-Wesley Professional, . Section 4.3.3.C: Discrete Fourier transforms, pg.305. The section begins by discussing the building blocks and concludes with a step-by-step description of the algorithm itself.

Convolutions

Suppose we are multiplying two numbers like 123 and 456 using long multiplication with base ''B'' digits, but without performing any carrying. The result might look something like this:

This sequence (4, 13, 28, 27, 18) is called the ''acyclic'' or ''linear convolution'' of the two original sequences (1,2,3) and (4,5,6). Once we have the acyclic convolution of two sequences, computing the product of the original numbers is easy: we just perform the carrying (for example, in the rightmost column, we keep the 8 and add the 1 to the column containing 27). In the example this yields the correct product 56088.

There are two other types of convolutions that will be useful. Suppose the input sequences have ''n'' elements (here 3). Then the acyclic convolution has ''n''+''n''−1 elements; if we take the rightmost ''n'' elements and add the leftmost ''n''−1 elements, this produces the cyclic convolution:

If we perform carrying on the cyclic convolution, the result is equivalent to the product of the inputs mod B^''n'' − 1. In the example, 10³ − 1 = 999, performing carrying on (28, 31, 31) yields 3141, and 3141 ≡ 56088 (mod 999).

Conversely, if we take the rightmost ''n'' elements and ''subtract'' the leftmost ''n''−1 elements, this produces the negacyclic convolution:

If we perform carrying on the negacyclic convolution, the result is equivalent to the product of the inputs mod B^''n'' + 1. In the example, 10³ + 1 = 1001, performing carrying on (28, 23, 5) yields 3035, and 3035 ≡ 56088 (mod 1001). The negacyclic convolution can contain negative numbers, which can be eliminated during carrying using borrowing, as is done in long subtraction.

Convolution theorem

Like other multiplication methods based on the fast Fourier transform, Schönhage–Strassen depends fundamentally on the convolution theorem, which provides an efficient way to compute the cyclic convolution of two sequences. It states that:

:The cyclic convolution of two vectors can be found by taking the discrete Fourier transform (DFT) of each of them, multiplying the resulting vectors element by element, and then taking the inverse discrete Fourier transform (IDFT).

Or in symbols:

:CyclicConvolution(''X, Y'') = IDFT(DFT(''X'') · DFT(''Y''))

If we compute the DFT and IDFT using a fast Fourier transform algorithm, and invoke our multiplication algorithm recursively to multiply the entries of the transformed vectors DFT(''X'') and DFT(''Y''), this yields an efficient algorithm for computing the cyclic convolution.

In this algorithm, it will be more useful to compute the ''negacyclic'' convolution; as it turns out, a slightly modified version of the convolution theorem (see discrete weighted transform) can enable this as well. Suppose the vectors X and Y have length ''n'', and ''a'' is a primitive root of unity of order 2''n'' (that is, ''a''^2''n'' = 1 and ''a'' to all smaller powers is not 1). Then we can define a third vector ''A'', called the ''weight vector'', as:

:''A'' = (''a''^''j''), 0 ≤ ''j'' < ''n''
:''A''⁻¹ = (''a''^−''j''), 0 ≤ ''j'' < ''n''

Now, we can state:

:NegacyclicConvolution(''X'', ''Y'') = ''A''⁻¹ · IDFT(DFT(''A'' · ''X'') · DFT(''A'' · ''Y''))

In other words, it's the same as before except that the inputs are first multiplied by ''A'', and the result is multiplied by ''A''⁻¹.

Choice of ring

The discrete Fourier transform is an abstract operation that can be performed in any algebraic ring; typically it's performed in the complex numbers, but actually performing complex arithmetic to sufficient precision to ensure accurate results for multiplication is slow and error-prone. Instead, we will use the approach of the number theoretic transform, which is to perform the transform in the integers mod ''N'' for some integer ''N''.

Just like there are primitive roots of unity of every order in the complex plane, given any order ''n'' we can choose a suitable N such that ''b'' is a primitive root of unity of order ''n'' in the integers mod ''N'' (in other words, ''b''^''n'' ≡ 1 (mod ''N''), and no smaller power of ''b'' is equivalent to 1 mod ''N'').

The algorithm will spend most of its time performing recursive multiplications of smaller numbers; with a naive algorithm, these occur in a number of places:

# Inside the fast Fourier transform algorithm, where the primitive root of unity ''b'' is repeatedly powered, squared, and multiplied by other values.
# When taking powers of the primitive root of unity ''a'' to form the weight vector A and when multiplying A or A⁻¹ by other vectors.
# When performing element-by-element multiplication of the transformed vectors.

The key insight to Schönhage–Strassen is to choose N, the modulus, to be equal to 2ⁿ + 1 for some integer ''n'' that is a multiple of the number of pieces ''P''=2^''p'' being transformed. This has a number of benefits in standard systems that represent large integers in binary form:

* Any value can be rapidly reduced modulo 2^''n'' + 1 using only shifts and adds, as explained in the next section.
* All roots of unity used by the transform can be written in the form 2^''k''; consequently we can multiply or divide any number by a root of unity using a shift, and power or square a root of unity by operating only on its exponent.
* The element-by-element recursive multiplications of the transformed vectors can be performed using a negacyclic convolution, which is faster than an acyclic convolution and already has "for free" the effect of reducing its result mod 2^''n'' + 1.

To make the recursive multiplications convenient, we will frame Schönhage–Strassen as being a specialized multiplication algorithm for computing not just the product of two numbers, but the product of two numbers mod 2^''n'' + 1 for some given ''n''. This is not a loss of generality, since one can always choose ''n'' large enough so that the product mod 2^''n'' + 1 is simply the product.

Shift optimizations

In the course of the algorithm, there are many cases in which multiplication or division by a power of two (including all roots of unity) can be profitably replaced by a small number of shifts and adds. This makes use of the observation that:

:(2^''n'')^k ≡ (−1)^k mod (2^''n'' + 1)

Note that a ''k''-digit number in base 2^''n'' written in positional notation can be expressed as (''d''_''k''−1,...,''d''₁,''d''₀). It represents the number $\scriptstyle \sum_^d_\cdot(2^)^$. Also note that for each ''d''_''i'' we have 0 ≤ ''d''_''i'' < 2^''n''.

This makes it simple to reduce a number represented in binary mod : take the rightmost (least significant) ''n'' bits, subtract the next ''n'' bits, add the next ''n'' bits, and so on until the bits are exhausted. If the resulting value is still not between 0 and 2^''n'', normalize it by adding or subtracting a multiple of the modulus . For example, if ''n''=3 (and so the modulus is 2³+1 = 9) and the number being reduced is 656, we have:

:656 = 1010010000₂ ≡ 000₂ − 010₂ + 010₂ − 1₂ = 0 − 2 + 2 − 1 = −1 ≡ 8 (mod 2³ + 1).