Schönhage-Strassen algorithm

The Schönhage-Strassen algorithm is an algorithm for multiplying two n -digit whole numbers . It was developed in 1971 by Arnold Schönhage and Volker Strassen . The algorithm is based on a very fast variant of the discrete, fast Fourier transformation as well as a clever change between the remainder class and cyclic arithmetic in finite number rings.

The Schönhage-Strassen algorithm terminates in (see Landau notation ) if the bit complexity on multi-volume Turing machines , i.e. the maximum runtime of the algorithm measured as required bit operations depending on the bit length of the input variables, is selected as the efficiency measure . This complexity represents an improvement both compared to the naive runtime algorithm known from school and compared to the Karatsuba algorithm developed in 1962 with a runtime of and its improved variant, the Toom Cook algorithm with runtime. ${\displaystyle O{\Big (}n\cdot \log(n)\cdot \log {\big (}\log(n){\big )}{\Big )}}$ ${\displaystyle n}$ ${\displaystyle O\left(n^{2}\right)}$ ${\displaystyle O\left(n^{\log _{2}(3)}\right)}$ ${\displaystyle O(n^{1+\varepsilon })}$

The Schönhage-Strassen algorithm was the most efficient known algorithm for the multiplication of large numbers from 1971 to 2007; In 2007 Martin Fürer published a further development of the algorithm with the even lower asymptotic complexity , where the iterated logarithm of n is. By optimizing Fürer's algorithm, David Harvey, Joris van der Hoeven and Grégoire Lecerf achieved a further improvement in the asymptotic running time . ${\displaystyle n\cdot \log(n)\cdot 2^{O(\log ^{*}(n))}}$ ${\displaystyle \log ^{*}(n)}$ ${\displaystyle O(n\cdot \log(n)\cdot 2^{3\log ^{*}(n)})}$

meaning

Until 2007, no more efficient algorithm could be found. In the general case, there is only the (trivial) linear running time as the lower bound, which the algorithm approaches with increasing number length. However, the researchers have found indications that the limit can never be undercut. Even with modern computers , this method of calculation is only more efficient than the Karatsuba algorithm for numbers with several thousand digits. This is probably less due to the overhead of the Schönhage-Strassen algorithm, but rather to the design optimization of computer processors that has been typical for decades , which gives preference to achieving fast floating point operations over arithmetic in finite residual class rings of whole numbers. ${\displaystyle O(n\cdot \log(n))}$

The Schönhage-Strassen algorithm is of central importance for the search for the algorithms with the best (time) complexity in computer algebra .

algorithm

Basic idea and terminology

The Schönhage – Strassen algorithm is based on the fast discrete Fourier transformation (DFT). This example shows the calculation of 1234 × 5678 = 7006652. The calculation takes place modulo 337. In order to improve the clarity, base 10 is used instead of base 2.

To multiply two whole numbers and , roughly the following scheme is used: ${\displaystyle a}$ ${\displaystyle b}$

1. Splitting the numbers (in binary representation) and into pieces of suitable length ${\displaystyle a}$ ${\displaystyle b}$
2. Fast discrete Fourier transformation (DFT) of the two piece sequences
3. Component-wise multiplication of the transformed pieces
4. Inverse transformation (inverse Fourier transformation) of the results
5. Combining the result pieces to form the result number

The small multiplications to be carried out in the middle step are in turn carried out in a recursive sense by the Schönhage-Strassen algorithm.

To understand why the result is the product of the numbers a and b, consider the polynomials

${\displaystyle A=\sum _{i=0}^{2^{n}-1}a_{i}\cdot X^{i}}$ and ${\displaystyle B=\sum _{i=0}^{2^{n}-1}b_{i}\cdot X^{i}}$

If you put in , you get the binary representation of the numbers a and b. To be calculated is for the product polynomial ${\displaystyle X=2}$ ${\displaystyle c=C(2)}$

${\displaystyle C=\sum _{l=0}^{2^{n+1}-1}\sum _{i=0}^{l}a_{i}b_{l-i}\cdot X^{l}}$

We determine the Fourier transform of the coefficient tuples of A and B:

${\displaystyle {\hat {a}}_{k}=\sum _{i=0}^{2^{n}-1}a_{i}\cdot w^{i\cdot k}}$ For ${\displaystyle k=0,\ldots ,2^{n}-1}$

${\displaystyle {\hat {b}}_{k}=\sum _{j=0}^{2^{n}-1}b_{j}\cdot w^{j\cdot k}}$ For ${\displaystyle k=0,\ldots ,2^{n}-1}$

In other words, you evaluate the two polynomials at the points . If these function values ​​are now multiplied, the corresponding function values ​​of the product polynomial result ${\displaystyle w^{k}}$

${\displaystyle C=A\cdot B}$.

To get the polynomial itself we need to undo the transformation: ${\displaystyle C}$

${\displaystyle {\hat {c}}_{l}=1/2^{n+1}\sum _{k=0}^{2^{n+1}-1}{\hat {a}}_{k}\cdot {\hat {b}}_{k}\cdot w^{-l\cdot k}}$ For ${\displaystyle l=0,\ldots ,2^{n+1}-1}$

${\displaystyle {\hat {c}}_{l}=1/2^{n+1}\sum _{i,j,k=0}^{2^{n+1}-1}a_{i}b_{j}\cdot w^{i\cdot k+j\cdot k-l\cdot k}}$ For ${\displaystyle l=0,\ldots ,2^{n+1}-1}$

${\displaystyle {\hat {c}}_{l}=1/2^{n+1}\sum _{i,j,k=0}^{2^{n+1}-1}a_{i}b_{j}\cdot w^{(i+j-l)\cdot k}}$ For ${\displaystyle l=0,\ldots ,2^{n+1}-1}$

After defining the roots of unity applies . This satisfies the following identity of geometric sums of roots of unity: ${\displaystyle w^{\,{2^{n+1}}}=1}$

${\displaystyle \displaystyle \sum _{k=0}^{{2^{n+1}}-1}w^{lk}={\begin{cases}2^{n+1}&l=0\\0&l\neq 0\end{cases}}}$ For ${\displaystyle l=0,...,2^{n+1}-1}$

because

${\displaystyle \displaystyle \sum _{k=0}^{2^{n+1}-1}x^{k}={\frac {x^{2^{n+1}}-1}{x-1}}}$ For ${\displaystyle x\neq 1}$

Thus:

${\displaystyle {\hat {c}}_{l}=\sum _{i+j=l}a_{i}b_{j}}$ For ${\displaystyle l=0,\ldots ,2^{n+1}-1}$

In the article Discrete Fourier Transformation , the mathematical basics of this transformation are detailed. Since the transformation results in sums with terms in each case , we still have a quadratic running time in a classic calculation of the terms ( e.g. using the Horner scheme ). These values ​​can be calculated faster using the fast Fourier transformation . This calculation is based on the following divide-and-conquer principle: ${\displaystyle 2^{n}}$ ${\displaystyle 2^{n}}$

${\displaystyle {\hat {a}}_{k}^{even}(w^{j})=\sum _{i=0}^{2^{n}-1}a_{j}\cdot w^{2j\cdot k}}$

${\displaystyle {\hat {a}}_{k}^{odd}(w^{j})=\sum _{i=0}^{2^{n}-1}a_{j}\cdot w^{(2j+1)\cdot k}}$

${\displaystyle {\hat {a}}_{k}(w^{j})=a_{k}^{even}(w^{2j})+w^{k}\cdot a_{k}^{odd}(w^{2j})}$

Partial solutions are put together using simple operations (addition and simple multiplication). This allows the transformations to be calculated in time . However, rounding the complex roots of units to a fixed length leads to calculation errors. In order to compensate for this, at least bits must be calculated for a resulting bit . This results in a total running time of . With the Schönhage-Strassen variant, we calculate in a remainder class ring instead and thus avoid the calculation errors of the complex numbers. ${\displaystyle O(N\cdot \log N)}$ ${\displaystyle O(\log N)}$ ${\displaystyle O(N\cdot (\log N)^{2})}$

Furthermore, the multiplication is not a pure convolution , but it can also lead to carryovers; after performing the FT and iFT, they must be treated appropriately.

The task of multiplying two whole numbers is now concretized as follows:

Let the two numbers to be multiplied be given in binary representation. Next is the maximum length (i.e. number of binary digits) of the two numbers. ${\displaystyle a,b\in \mathbb {Z} }$ ${\displaystyle N}$

After appropriate treatment of the signs of the two numbers as well as the trivial special cases and (which can be done with linear effort ) one can assume that numbers are natural. The Schönhage-Strassen algorithm solves this task in . ${\displaystyle a=0}$ ${\displaystyle b=0}$ ${\displaystyle O(N)}$ ${\displaystyle a,b\in \mathbb {N} }$ ${\displaystyle O(N\cdot \log N\cdot \log \log N)}$

Theoretical preparation

Super fast DFT

The super-fast DFT mentioned above, which is the core of the algorithm, needs to be explained in more detail, as it is used here very specifically.

Let it be a commutative unitary ring . In let the element be a unit ; continue to be a th root of unity (thus ) which fulfills equality . Then the computation of the discrete Fourier transform (DFT) in the product space (this is a short notation for ; the term vector space is only used here for the case that there is a body ) can be carried out in a fast variant (as FFT) as follows: ${\displaystyle R}$ ${\displaystyle R}$ ${\displaystyle 2}$ ${\displaystyle w\in R}$ ${\displaystyle 2^{n}}$ ${\displaystyle w^{2^{n}}=1}$ ${\displaystyle w^{2^{n-1}}=-1}$ ${\displaystyle R^{2^{n}}}$ ${\displaystyle R^{(2^{n})}}$ ${\displaystyle R}$

Calculate for the transformed with ${\displaystyle a=(a_{0},\ldots ,a_{2^{n}-1})\in R^{2^{n}}}$ ${\displaystyle {\hat {a}}\in R^{2^{n}}}$

${\displaystyle {\hat {a}}_{k}=\sum _{j=0}^{2^{n}-1}a_{j}\cdot w^{j\cdot k}}$for . ${\displaystyle k=0,\ldots ,2^{n}-1}$

By writing down the indices and in binary representation, doing this in reverse order for the number , the transform can be optimally calculated as follows: ${\displaystyle k=\sum _{\nu =0}^{n-1}k_{\nu }\cdot 2^{\nu }}$ ${\displaystyle j=\sum _{\nu =0}^{n-1}j_{\nu }\cdot 2^{n-1-\nu }}$ ${\displaystyle j}$ ${\displaystyle {\hat {a}}}$

Be there

${\displaystyle A_{0}(j_{0},\ldots j_{n-1})=a_{j}}$ For ${\displaystyle j=0,\ldots ,2^{n}-1}$

and

${\displaystyle A_{r+1}(k_{0},\ldots ,k_{r},j_{r+1},\ldots ,j_{n-1})=}$

${\displaystyle =A_{r}(k_{0},\ldots ,k_{r-1},0,j_{r+1},\ldots ,j_{n-1})+}$

${\displaystyle A_{r}(k_{0},\ldots ,k_{r-1},1,j_{r+1},\ldots ,j_{n-1})\cdot w^{2^{n-1-r}\cdot (k_{r}2^{r}+\dots +k_{0}2^{0})}}$

${\displaystyle =\sum _{j_{r}=0}^{1}A_{r}(k_{0},\ldots ,k_{r-1},j_{r},\ldots ,j_{n-1})\cdot w^{j_{r}2^{n-1-r}\cdot (k_{r}2^{r}+\dots +k_{0}2^{0})}}$.

The closed representation for these intermediate terms is

${\displaystyle A_{r+1}(k_{0},\ldots ,k_{r},j_{r+1},\ldots ,j_{n-1})=}$

${\displaystyle =\sum _{j_{r}=0}^{1}\ldots \sum _{j_{0}=0}^{1}a_{j}\cdot w^{(j_{0}2^{n-1}+\ldots +j_{r}2^{n-1-r})\cdot (k_{r}2^{r}+\dots +k_{0}2^{0})}}$.

(To recalculate this representation, note ). ${\displaystyle w^{(j_{0}2^{n-1}+\ldots +j_{r-1}2^{n-r})\cdot k_{r}2^{r}}=1}$

This recursion provides the desired Fourier coefficients . ${\displaystyle {\hat {a}}_{k}=A_{n}(k_{0},\ldots ,k_{n-1})}$

Due to the property , we can transform the recursion step to be somewhat more computationally friendly ${\displaystyle w^{2^{n-1}}=-1}$

${\displaystyle A_{r+1}(k_{0},\ldots ,k_{r-1},0,j_{r+1},\ldots ,j_{n-1})=}$

${\displaystyle =A_{r}(k_{0},\ldots ,k_{r-1},0,j_{r+1},\ldots ,j_{n-1})+A_{r}(k_{0},\ldots ,k_{r-1},1,j_{r+1},\ldots ,j_{n-1})\cdot w^{x}}$

and

${\displaystyle A_{r+1}(k_{0},\ldots ,k_{r-1},1,j_{r+1},\ldots ,j_{n-1})=}$

${\displaystyle =A_{r}(k_{0},\ldots ,k_{r-1},0,j_{r+1},\ldots ,j_{n-1})-A_{r}(k_{0},\ldots ,k_{r-1},1,j_{r+1},\ldots ,j_{n-1})\cdot w^{x}}$

with the same exponent . ${\displaystyle x=2^{n-1-r}\cdot (k_{r-1}2^{r-1}+\ldots +k_{0}2^{0})}$

The inverse transformation, i.e. the inverse FFT, succeeds because we have assumed that it is invertible in the ring : ${\displaystyle 2}$ ${\displaystyle R}$

${\displaystyle A_{r}(k_{0},\ldots ,k_{r-1},0,j_{r+1},\ldots ,j_{n-1})=}$

${\displaystyle =2^{-1}{\big (}A_{r+1}(k_{0},\ldots ,k_{r-1},0,j_{r+1},\ldots ,j_{n-1})+A_{r+1}(k_{0},\ldots ,k_{r-1},1,j_{r+1},\ldots ,j_{n-1})}$

such as

${\displaystyle A_{r}(k_{0},\ldots ,k_{r-1},1,j_{r+1},\ldots ,j_{n-1})=}$

${\displaystyle =2^{-1}w^{-x}{\big (}A_{r+1}(k_{0},\ldots ,k_{r-1},0,j_{r+1},\ldots ,j_{n-1})-A_{r+1}(k_{0},\ldots ,k_{r-1},1,j_{r+1},\ldots ,j_{n-1})}$,

where again is. ${\displaystyle x=2^{n-1-r}\cdot (k_{r-1}2^{r-1}+\ldots +k_{0}2^{0})}$

When used in the Schönhage-Strassen algorithm, only a halved FFT is actually required; This means the following: Let us start with the calculation in the first step of the recursion

${\displaystyle A_{1}(1,j_{1},\ldots ,j_{n-1})=A_{0}(0,j_{1},\ldots ,j_{n-1})-A_{0}(1,j_{1},\ldots ,j_{n-1})=a_{j}-a_{j+2^{n-1}}}$

only for and if we also restrict the further steps of the recursion to , we calculate all of them for odd values . Conversely, if one only wants to recover the differences from the original ones for odd (these are pieces) , then the recursion halved in the reverse direction is sufficient. ${\displaystyle k_{0}=1}$ ${\displaystyle k_{0}=1}$ ${\displaystyle {\hat {a}}_{k}}$ ${\displaystyle k}$ ${\displaystyle {\hat {a}}_{k}}$ ${\displaystyle k}$ ${\displaystyle 2^{n-1}}$ ${\displaystyle a_{j}-a_{j+2^{n-1}}}$ ${\displaystyle a_{j}}$

In the Schönhage-Strassen algorithm, the described fast Fourier transformation is required for finite number rings with Fermat numbers . ${\displaystyle \mathbb {Z} _{F_{n}}}$ ${\displaystyle F_{n}=2^{(2^{n})}+1}$

Note on the notation: For the remainder class ring we use the shorter notation , which could only lead to confusion in the context of the p-adic numbers . ${\displaystyle \mathbb {Z} /k\mathbb {Z} }$ ${\displaystyle \mathbb {Z} _{k}}$

The number (or, depending on the context, a suitable power of 2) will be used as the root of unity in the ring . The multiplications to be carried out with the FFT algorithm are then of the form ; however, they cannot be carried out as pure shift operations, since the reduction of a larger intermediate result has to be postponed modulo . This is where one of the brilliant ideas from Schönhage and Strassen comes into play : They embed the ring (equipped with the residual class arithmetic ) in a larger ring equipped with the cyclic arithmetic . This overring has a power of 2 as order, so that the corresponding multiplication can actually be carried out in it as a pure shift operation. This trick can be summarized in a nice structural theorem about residual class and cyclic arithmetic in finite number rings. ${\displaystyle \mathbb {Z} _{F_{n}}}$ ${\displaystyle 2}$ ${\displaystyle x\cdot 2^{k}}$ ${\displaystyle F_{n}}$

Structural theorem about cyclic arithmetic

The structural theorem about cyclic arithmetic can be formally expressed as follows:

For a power of two with a natural number, the following applies ${\displaystyle D=2^{n}}$ ${\displaystyle n\in \mathbb {N} }$

${\displaystyle (\mathbb {Z} _{D+1},+,\cdot )\cong (\mathbb {Z} _{D^{2}},\oplus ,\otimes )/(D+1)\cdot \mathbb {Z} }$.

Here referred by the representatives displayable residue classes modulo equipped with the residual arithmetic, d. H. with addition and multiplication modulo . The numbers occurring in this remainder class ring can be represented with binary digits. ${\displaystyle (\mathbb {Z} _{D+1},+,\cdot )}$ ${\displaystyle 0,\ldots ,D}$ ${\displaystyle D+1}$ ${\displaystyle D+1}$ ${\displaystyle n+1}$

The structure on the right-hand side denotes the remainder classes modulo of the number , which, however, are not equipped with the remainder class arithmetic, but instead with the cyclic arithmetic. In the case of intermediate results that are too large, carry-overs are canceled and added to the final result. In the binary digit representation, this corresponds to a shift of the excess binary digits (right-aligned at the lowest digit positions) with subsequent addition. For example, adding with does not give the value , but the value . From the number structure obtained in this way with cyclic arithmetic, the factor ring is now formed modulo . So the end results are reduced modulo . ${\displaystyle (\mathbb {Z} _{D^{2}},\oplus ,\otimes )}$ ${\displaystyle D^{2}}$ ${\displaystyle a+1}$ ${\displaystyle a=D^{2}-1}$ ${\displaystyle D^{2}\equiv 0}$ ${\displaystyle D^{2}+1\equiv 1}$ ${\displaystyle D+1}$ ${\displaystyle D+1}$

This structural theorem thus means the following: The modulo calculation in can also be replaced by the cyclic calculation in the larger number space with subsequent reduction modulo . ${\displaystyle \mathbb {Z} _{D+1}}$ ${\displaystyle \mathbb {Z} _{D^{2}}}$ ${\displaystyle D+1}$

The decisive factor for the success of the embedding presented in this structural theorem is the property that the largest representable number in the cyclic number space (here this is the number ) represents the number from the remainder class ring . The condition is necessary for this. On the other hand, for cyclic arithmetic to be meaningfully defined, it must be a power of two. Together it results that the optimal choice for the size of the cyclic embedding space is. ${\displaystyle F}$ ${\displaystyle D^{2}-1}$ ${\displaystyle 0}$ ${\displaystyle \mathbb {Z} _{D+1}}$ ${\displaystyle (D+1)|F}$ ${\displaystyle F+1}$ ${\displaystyle F=D^{2}-1}$

The classic residual class ring , on the other hand, would not be suitable for embedding, because in this ring the following applies : H. the number in this ring is a zero divisor . ${\displaystyle (\mathbb {Z} _{D^{2}},+,\cdot )}$ ${\displaystyle 2^{2n}=D^{2}=0}$ ${\displaystyle 2}$

execution

If we have the numbers to be multiplied with binary digits, we carry out different recursion steps, depending on whether it is even or odd, in order to log the number of digits in a single step: ${\displaystyle a,b}$ ${\displaystyle N\leq 2^{m}}$ ${\displaystyle m}$

Recursion step for odd m

We carry out this step of tracing back from to with the complexity . ${\displaystyle m=2n-1}$ ${\displaystyle n}$ ${\displaystyle O\left(2^{n}\psi (2^{n})+n\cdot 2^{2n}\right)}$

It is to be multiplied by and the Fermat number . In this step we will carry out the return to the Fermat number . ${\displaystyle a,b\in \mathbb {Z} _{F_{m}}}$ ${\displaystyle m=2n-1}$ ${\displaystyle F_{m}=2^{2^{m}}+1}$ ${\displaystyle F_{n}=2^{2^{n}}+1}$

We use the abbreviations for the powers of two belonging to the two Fermat numbers

${\displaystyle D=F_{n}-1=2^{2^{n}}}$

and

${\displaystyle E=F_{m}-1=2^{2^{2n-1}}=D^{2^{n-1}}}$

a. Halving the number of digits will become our denomination size; H. we develop and according to potencies of : ${\displaystyle D}$ ${\displaystyle a}$ ${\displaystyle b}$ ${\displaystyle {\sqrt {D}}}$

${\displaystyle a=\sum _{i=0}^{2^{n}}a_{i}\cdot {\sqrt {D}}^{i}\quad }$and , ${\displaystyle \quad b=\sum _{i=0}^{2^{n}}b_{i}\cdot {\sqrt {D}}^{i}}$

where applies to the individual pieces . In binary representation, this breakdown corresponds to a simple grouping of the bit sequences into pieces of bits length . ${\displaystyle 0\leq a_{i},b_{i}<{\sqrt {D}}}$ ${\displaystyle 2^{n-1}}$

A small weakness of the algorithm (which, however, does not detract from the complexity limit reached) is now revealed. In order to be able to use the super-fast DFT on the piece sequences and , these must be padded with zeros to the next power of two length; the number representation is thus artificially extended to ${\displaystyle (a_{0},\ldots ,a_{2^{n}})}$ ${\displaystyle (b_{0},\ldots ,b_{2^{n}})}$

${\displaystyle a=\sum _{i=0}^{2^{n+1}-1}a_{i}\cdot {\sqrt {D}}^{i}}$and . ${\displaystyle b=\sum _{j=0}^{2^{n+1}-1}b_{j}\cdot {\sqrt {D}}^{j}}$

By virtue of the structural theorem on cyclic arithmetic mentioned above, we now switch from the remainder class ring to the quotient space with cyclic arithmetic. In this space is calculated for the multiplication problem ${\displaystyle \mathbb {Z} _{E+1}}$ ${\displaystyle (\mathbb {Z} _{E^{2}},\oplus ,\otimes )/(E+1)\cdot \mathbb {Z} }$

${\displaystyle a\cdot b=\left(\sum _{i=0}^{2^{n+1}-1}a_{i}{\sqrt {D}}^{i}\right)\cdot \left(\sum _{j=0}^{2^{n+1}-1}b_{j}{\sqrt {D}}^{j}\right)}$

${\displaystyle =\sum _{k=0}^{2^{n+2}-2}\left(\sum _{i,j=0 \atop i+j=k}^{2^{n+1}-1}a_{i}\cdot b_{j}\right){\sqrt {D}}^{k}}$

${\displaystyle =\sum _{k=0}^{2^{n+1}-1}\left(\sum _{i,j=0 \atop i+j=k}^{2^{n+1}-1}a_{i}\cdot b_{j}\right){\sqrt {D}}^{k}+\sum _{k=0}^{2^{n+1}-1}\left(\sum _{i,j=0 \atop i+j=2^{n+1}+k}^{2^{n+1}-1}a_{i}\cdot b_{j}\right){\sqrt {D}}^{2^{n+1}+k}}$

${\displaystyle =\sum _{k=0}^{2^{n+1}-1}\left(\sum _{i,j=0 \atop i+j\equiv k\mod 2^{n+1}}^{2^{n+1}-1}a_{i}\cdot b_{j}\right){\sqrt {D}}^{k}}$,

where in the last step we used the property in this cyclic number space. ${\displaystyle {\sqrt {D}}^{2^{n+1}}=D^{2^{n}}=E^{2}=1}$

In summary, the multiplication gets the form

${\displaystyle a\cdot b=\sum _{k=0}^{2^{n+1}-1}c_{k}{\sqrt {D}}^{k}}$

with the result coefficients

${\displaystyle c_{k}=\sum _{i,j=0 \atop i+j\equiv k\mod 2^{n+1}}^{2^{n+1}-1}a_{i}\cdot b_{j}}$.

We can estimate upwards. ${\displaystyle c_{k}<2^{n+1}D}$

The sum formula is now described so that we can limit the FFT to be applied to a halved FFT.

It applies , therefore is ${\displaystyle c_{k+2^{n}}{\sqrt {D}}^{k+2^{n}}=c_{k+2^{n}}{\sqrt {D}}^{k}E=-c_{k+2^{n}}{\sqrt {D}}^{k}}$

${\displaystyle a\cdot b=\sum _{k=0}^{2^{n}-1}(c_{k}-c_{k+2^{n}}){\sqrt {D}}^{k}}$

with in . With suitable addition we can shift the range of values ​​into the positive, namely it is , and with the definition ${\displaystyle -2^{n+1}D<c_{k}-c_{k+2^{n}}<2^{n+1}D}$ ${\displaystyle \mathbb {Z} _{E+1}}$ ${\displaystyle 0<c_{k}-c_{k+2^{n}}+2^{n+1}D<2^{n+2}D}$

${\displaystyle z_{k}={\Big \lbrace }{c_{k}-c_{k+2^{n}}+2^{n+1}D,\quad 0\leq k<2^{n} \atop 2^{n+1}D,\quad \quad 2^{n}\leq k<2^{n+1}}}$

applies

${\displaystyle a\cdot b=\sum _{k=0}^{2^{n+1}-1}z_{k}{\sqrt {D}}^{k}}$.

The estimate applies to the non-trivial ones (indices to ) . Since the two numbers and are relatively prime, it is sufficient to calculate the remainders and to determine the . ${\displaystyle z_{k}}$ ${\displaystyle 0}$ ${\displaystyle 2^{n}-1}$ ${\displaystyle 0<z_{k}<2^{n+2}D<2^{n+2}F_{n}}$ ${\displaystyle 2^{n+2}}$ ${\displaystyle F_{n}}$ ${\displaystyle z_{k}}$ ${\displaystyle z_{k}\mod 2^{n+2}}$ ${\displaystyle z_{k}\mod F_{n}}$

If you have determined the remainder and , you can calculate in complexity as follows: Calculate first and then . ${\displaystyle z_{k}=\xi \mod F_{n}}$ ${\displaystyle z_{k}=\eta \mod 2^{n+2}}$ ${\displaystyle O(2^{n})}$ ${\displaystyle \delta =\eta -\xi \mod 2^{n+2}}$ ${\displaystyle z_{k}=\xi +\delta \cdot (D+1)=\xi +\delta \cdot F_{n}}$

Determination of the remainders modulo 2 ^{n + 2}

Here we use a trick that is very typical for computer algebra: We put the piece sequences and insert sufficiently long zero sequences with safety margins so that, after the product has been formed, the individual results are also lined up in pieces without overlapping. So be and in . We are now forming ${\displaystyle a_{i}}$ ${\displaystyle b_{i}}$ ${\displaystyle \alpha _{j}=a_{j}\mod 2^{n+2}}$ ${\displaystyle \beta _{j}=b_{j}\mod 2^{n+2}}$ ${\displaystyle \mathbb {Z} _{2^{n+2}}}$

${\displaystyle u=\sum _{k=0}^{2^{n+1}-1}\alpha _{k}2^{k(3n+5)}}$ and ${\displaystyle v=\sum _{k=0}^{2^{n+1}-1}\beta _{k}2^{k(3n+5)}}$

and have with us . The product then contains the sums in disjoint pieces of the bit length ${\displaystyle 0\leq u,v<2^{2^{n+1}(3n+5)}}$ ${\displaystyle u\cdot v}$ ${\displaystyle 3n+5}$

${\displaystyle \gamma _{k}=\sum _{r,s=0 \atop r+s=k}^{2^{n+1}-1}\alpha _{r}\cdot \beta _{s}}$

with because it is . For the terms of our original multiplication problem we see ${\displaystyle 0\leq k<2^{n+2}}$ ${\displaystyle 0\leq \gamma _{k}<2^{3n+5}}$ ${\displaystyle c_{k}}$ ${\displaystyle a\cdot b}$

${\displaystyle c_{k}=\gamma _{k}+\gamma _{k+2^{n+1}}\mod 2^{n+2}}$.

For the remnants to be determined we get ${\displaystyle \eta _{k}=z_{k}\mod 2^{n+2}}$

${\displaystyle \eta _{k}=\gamma _{k}-\gamma _{k+2^{n}}+\gamma _{k+2\cdot 2^{n}}-\gamma _{k+3\cdot 2^{n}}}$in . ${\displaystyle \mathbb {Z} _{2^{n+2}}}$

The complexity involved in forming all as well as extracting the is ; multiplication costs , so this is total . ${\displaystyle \alpha _{j},\beta _{j}}$ ${\displaystyle \eta _{k}}$ ${\displaystyle O(2^{2n})}$ ${\displaystyle u\cdot v}$ ${\displaystyle \psi (2^{n+1}(3n+5))}$ ${\displaystyle O(2^{2n})}$

Determination of the remainders modulo (D + 1)

This is where the DFT comes in. We subject the vectors and with the DFT in with and the number as the -th root of unity. Since we only need the differences , the halved DFT is sufficient: ${\displaystyle (a_{0},\ldots a_{2^{n+1}-1})}$ ${\displaystyle (b_{0},\ldots b_{2^{n+1}-1})}$ ${\displaystyle 0\leq a_{k},b_{k}<{\sqrt {D}}}$ ${\displaystyle R^{2^{n+1}}}$ ${\displaystyle R=\mathbb {Z} _{D+1}}$ ${\displaystyle 2}$ ${\displaystyle 2^{n+1}}$ ${\displaystyle c_{k}-c_{k+2^{n}}}$

• DFT to determine the and only for the odd with ${\displaystyle {\hat {a}}_{k}}$ ${\displaystyle {\hat {b}}_{k}}$ ${\displaystyle k}$ ${\displaystyle 0\leq k<2^{n+1}}$
• ${\displaystyle 2^{n}}$Multiplications for all odd ${\displaystyle {\hat {c}}_{k}={\hat {a}}_{k}\cdot {\hat {b}}_{k}}$ ${\displaystyle k}$
• Inverse DFT to get all the differences from the odd ${\displaystyle c_{k}-c_{k+2^{n}}}$ ${\displaystyle {\hat {c}}_{k}}$ ${\displaystyle k}$

The complexity outlay for this consists of steps of the individual outlay for the DFT (ie total ); Added to this are the addition of and the reductions modulo for the production of what can be mastered in. ${\displaystyle O(n2^{n})}$ ${\displaystyle O(2^{n})}$ ${\displaystyle O(n2^{2n})}$ ${\displaystyle 2^{n+1}D}$ ${\displaystyle (D+1)}$ ${\displaystyle z_{k}\mod (D+1)}$ ${\displaystyle O(2^{2n})}$

Recursion step for even m

The complexity is also achieved for this step of tracing back from to . ${\displaystyle m=2n-2}$ ${\displaystyle n}$ ${\displaystyle O\left(2^{n}\psi (2^{n})+n\cdot 2^{2n}\right)}$

It is to be multiplied by and the Fermat number . In this step, too, we will return to the Fermat number . ${\displaystyle a,b\in \mathbb {Z} _{F_{m}}}$ ${\displaystyle m=2n-2}$ ${\displaystyle F_{m}=2^{2^{m}}+1}$ ${\displaystyle F_{n}=2^{2^{n}}+1}$

We use the same abbreviations for the powers of two belonging to the two Fermat numbers

${\displaystyle D=F_{n}-1=2^{2^{n}}}$

and

${\displaystyle E=F_{m}-1=2^{2^{2n-2}}=D^{2^{n-2}}}$

a. Again, the halved number of digits will become our denomination size; H. we develop and according to potencies of : ${\displaystyle D}$ ${\displaystyle a}$ ${\displaystyle b}$ ${\displaystyle {\sqrt {D}}}$

${\displaystyle a=\sum _{i=0}^{2^{n-1}}a_{i}\cdot {\sqrt {D}}^{i}\quad }$and , ${\displaystyle \quad b=\sum _{i=0}^{2^{n-1}}b_{i}\cdot {\sqrt {D}}^{i}}$

where applies to the individual pieces . ${\displaystyle 0\leq a_{i},b_{i}<{\sqrt {D}}}$

As above, we extend the number representation to a power of two

${\displaystyle a=\sum _{i=0}^{2^{n}-1}a_{i}\cdot {\sqrt {D}}^{i}}$

and analog for . ${\displaystyle b}$

With the help of the structural theorem for cyclic arithmetic, we now switch from the remainder class ring to the quotient space with cyclic arithmetic. ${\displaystyle \mathbb {Z} _{E+1}}$ ${\displaystyle (\mathbb {Z} _{E^{2}},\oplus ,\otimes )/(E+1)\cdot \mathbb {Z} }$

We can do that again

${\displaystyle a\cdot b=\sum _{k=0}^{2^{n}-1}c_{k}{\sqrt {D}}^{k}}$

with the result coefficients

${\displaystyle c_{k}=\sum _{r,s=0 \atop r+s\equiv k\mod 2^{n}}^{2^{n}-1}a_{r}\cdot b_{s}}$

represent. We can estimate upwards. ${\displaystyle c_{k}<2^{n}D}$

For we can again ${\displaystyle {\sqrt {D}}^{2^{n-1}}=E^{2}=1}$

${\displaystyle 0<c_{k}-c_{k+2^{n-1}}+2^{n}D<2^{n+1}D}$

infer, and with

${\displaystyle z_{k}={\Big \lbrace }{c_{k}-c_{k+2^{n-1}}+2^{n}D,\quad 0\leq k<2^{n-1} \atop 2^{n}D,\quad \quad 2^{n-1}\leq k<2^{n}}}$

applies

${\displaystyle a\cdot b=\sum _{k=0}^{2^{n}-1}z_{k}{\sqrt {D}}^{k}}$

with . The estimate applies to the non-trivial ones (indices to ) . Because of the coprime nature of the two numbers and it is again sufficient to determine the , the remainders and to calculate. ${\displaystyle 0<z_{k}<2^{n+1}D}$ ${\displaystyle z_{k}}$ ${\displaystyle 0}$ ${\displaystyle 2^{n-1}-1}$ ${\displaystyle 0<z_{k}<2^{n+1}D<2^{n+1}F_{n}}$ ${\displaystyle 2^{n+1}}$ ${\displaystyle F_{n}}$ ${\displaystyle z_{k}}$ ${\displaystyle z_{k}\mod 2^{n+2}}$ ${\displaystyle z_{k}\mod F_{n}}$

Determination of the remainder modulo 2 ^{n + 1}

We use the trick of inserting safety margins again: Let it be and in . We educate ${\displaystyle \alpha _{j}=a_{j}\mod 2^{n+1}}$ ${\displaystyle \beta _{j}=b_{j}\mod 2^{n+1}}$ ${\displaystyle \mathbb {Z} _{2^{n+1}}}$

${\displaystyle u=\sum _{k=0}^{2^{n}-1}\alpha _{k}2^{k(3n+2)}}$ and ${\displaystyle v=\sum _{k=0}^{2^{n}-1}\beta _{k}2^{k(3n+2)}}$

and have with us . The product then contains the sums in disjoint pieces of the bit length ${\displaystyle 0\leq u,v<2^{2^{n}(3n+2)}}$ ${\displaystyle u\cdot v}$ ${\displaystyle 3n+2}$

${\displaystyle \gamma _{k}=\sum _{r,s=0 \atop r+s=k}^{2^{n}-1}\alpha _{r}\cdot \beta _{s}}$

with . For our original multiplication problem we are looking for ${\displaystyle 0\leq k<2^{n+1}}$ ${\displaystyle c_{k}}$ ${\displaystyle a\cdot b}$

${\displaystyle c_{k}=\gamma _{k}+\gamma _{k+2^{n}}\mod 2^{n+1}}$.

For the remnants to be determined we get ${\displaystyle \eta _{k}=z_{k}\mod 2^{n+1}}$

${\displaystyle \eta _{k}=\gamma _{k}-\gamma _{k+2^{n-1}}+\gamma _{k+2\cdot 2^{n-1}}-\gamma _{k+3\cdot 2^{n-1}}}$in . ${\displaystyle \mathbb {Z} _{2^{n+1}}}$

Determination of the remainders modulo (D + 1)

With we subject the vectors again and with the DFT in , this time choosing the number as the -th root of unity. Since we only need the differences , the halved DFT is sufficient here: ${\displaystyle R=\mathbb {Z} _{D+1}}$ ${\displaystyle (a_{0},\ldots a_{2^{n}-1})}$ ${\displaystyle (b_{0},\ldots b_{2^{n}-1})}$ ${\displaystyle 0\leq a_{k},b_{k}<{\sqrt {D}}}$ ${\displaystyle R^{2^{n}}}$ ${\displaystyle 4}$ ${\displaystyle 2^{n}}$ ${\displaystyle c_{k}-c_{k+2^{n-1}}}$

• DFT to determine the and only for the odd with ${\displaystyle {\hat {a}}_{k}}$ ${\displaystyle {\hat {b}}_{k}}$ ${\displaystyle k}$ ${\displaystyle 0\leq k<2^{n}}$
• ${\displaystyle 2^{n-1}}$Multiplications for all odd ${\displaystyle {\hat {c}}_{k}={\hat {a}}_{k}\cdot {\hat {b}}_{k}}$ ${\displaystyle k}$
• Inverse DFT to get all the differences from the odd ${\displaystyle c_{k}-c_{k+2^{n-1}}}$ ${\displaystyle {\hat {c}}_{k}}$ ${\displaystyle k}$

Summary

Starting with and with digit length , the recursion shown achieves a complexity of . ${\displaystyle a}$ ${\displaystyle b}$ ${\displaystyle n}$ ${\displaystyle O(n\cdot \log(n)\cdot \log(\log(n)))}$

Modified form

Zimmermann and Brent describe a variant of the algorithm in which the runtime (depending on the length of the input) does not make jumps, but rather runs more steadily. This is achieved in that the DFT vectors are not formed from -digit binary numbers, but numbers of the appropriate length. This means that the length of the vectors to be transformed does not have to be a power of two. ${\displaystyle 2^{n}}$

literature

• Arnold Schönhage : Asymptotically fast algorithms for the numerical multiplication and division of polynomials with complex coefficients . In: J. Calmet. (Ed.): EUROCAM '82: European Computer Algebra Conference (Marseille, France, April 1982) . Lect. Notes Comp. Sci. 144. Springer, 1982.
• Donald E. Knuth : The Art of Computer Programming , Vol. 2: Seminumerical Algorithms . 3. Edition. Addison-Wesley, 1998, ISBN 0-201-89684-2 .
• Chee Yap, Chen Li: QuickMul: Practical FFT-based Integer Multiplication , 2000. Simplification of the Schönhage-Strassen algorithm for practical applications.
• Michael T. Goodrich, Roberto Tamassia: Algorithm Design Foundations, Analysis, and Internet Examples (PDF; 308 kB). An introduction to FFT with a Java implementation of the QuickMul algorithm.
• Daniel J. Bernstein : Multidigit multiplication for mathematicians . August 11, 2001 ( yp.to - summary of different techniques for polynomial and long number multiplication ).

Web links

• World record calculation method comes in late honors . University of Bonn, press release, December 21, 2004

Individual evidence

1. ^ Arnold Schönhage, Volker Strassen : Fast multiplication of large numbers . In: Computing , 7, 1971, pp. 281-292, Springer Verlag
2. ↑ Martin Fürer: Faster integer multiplication . ( Memento of the original from April 25, 2013 in the Internet Archive ; PDF; 295 kB) Info: The archive link was inserted automatically and has not yet been checked. Please check the original and archive link according to the instructions and then remove this notice. STOC 2007 Proceedings, pp. 57-66. 
3. ^ David Harvey, Joris van der Hoeven, Grégoire Lecerf: Even faster integer multiplication . 2014, arxiv : 1407.3360
4. ↑ loria.fr (PDF) p. 56
5. ↑ loria.fr (PDF)