11-05-2012, 11:10 AM
Schonhage–Strassen algorithm
Strassen algorithm.docx (Size: 61.58 KB / Downloads: 40)
HISTORY
The Schonhage–Strassen algorithm is an asymptotically fast multiplication for large integers. It was developed by Arnold Schönhage and Volker Strassen in 1971. The run-time bit complexity is, in Big O notation, O(N log N log log N), while the arithmetic complexity is O(N log N). The algorithm uses recursive Fast Fourier transforms in rings with 22n + 1 element, a specific type of number theoretic transform.
The Schonhage–Strassen algorithm was the asymptotically fastest multiplication method known from 1971 until 2007, when a new method, Furer's algorithm, was announced with lower asymptotic complexity; however, Furer's algorithm currently only achieves an advantage for astronomically large values and is not used in practice.
In practice the Schonhage–Strassen algorithm starts to outperform older methods such as Karatsuba and Toom–Cook multiplication for numbers beyond 2215 to 2217 (10,000 to 40,000 decimal digits).
Details
This section explains in detail how Schönhage–Strassen is implemented. It is based primarily on an overview of the method by Crandall and Pomerance in their Prime Numbers: A Computational Perspective. Another source for detailed information is Knuth's The Art of Computer Programming.
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:
1 2 3
× 4 5 6
________________________________________
6 12 18
5 10 15
4 8 12
________________________________________
4 13 28 27 18
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 you have the acyclic convolution of two sequences, computing the product of the original numbers is easy: you just perform the carrying (for example, in the rightmost column, you'd 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:
28 27 18
+ 4 13
________________________________________
28 31 31
If we perform carrying on the cyclic convolution, the result is equivalent to the product of the inputs mod Bn − 1. In the example, 103 − 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:
28 27 18
− 4 13
________________________________________
28 23 5
If we perform carrying on the negacyclic convolution, the result is equivalent to the product of the inputs mod Bn + 1. In the example, 103 + 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 can enable this as well. Suppose the vectors X and Y have length n, and a is a primitive root of unity of order 2n (that is, a2n = 1 and a to all smaller powers is not 1). Then we can define a third vector A, called the weight vector, as:
A = (aj), 0 ≤ j < n
A−1 = (a−j), 0 ≤ j < n
Now, we can state:
NegacyclicConvolution(X, Y) = A−1 • 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−1.