This articleneeds additional citations forverification. Please helpimprove this article byadding citations to reliable sources. Unsourced material may be challenged and removed. Find sources: "Schönhage–Strassen algorithm" – news ·newspapers ·books ·scholar ·JSTOR(October 2024) (Learn how and when to remove this message)

TheSchönhage–Strassen algorithm is an asymptotically fastmultiplication algorithm for largeintegers, published byArnold Schönhage andVolker Strassen in 1971.^[1] It works by recursively applyingfast Fourier transform (FFT) overthe integers modulo${\displaystyle 2^n+1}$. The run-timebit complexity to multiply twon-digit numbers using the algorithm is${\displaystyle O(n\cdot \log n\cdot \log \log n)}$ inbigO notation.

The Schönhage–Strassen algorithm was theasymptotically fastest multiplication method known from 1971 until 2007. It is asymptotically faster than older methods such asKaratsuba andToom–Cook multiplication, and starts to outperform them in practice for numbers beyond about 10,000 to 100,000 decimal digits.^[2] In 2007, Martin Fürer publishedan algorithm with faster asymptotic complexity.^[3] In 2019, David Harvey andJoris van der Hoeven demonstrated that multi-digit multiplication has theoretical${\displaystyle O(n\log n)}$ complexity; however, their algorithm has constant factors which make it impossibly slow for any conceivable practical problem (seegalactic algorithm).^[4]

Applications of the Schönhage–Strassen algorithm include large computations done for their own sake such as theGreat Internet Mersenne Prime Search andapproximations ofπ, as well as practical applications such asLenstra elliptic curve factorization viaKronecker substitution, which reduces polynomial multiplication to integer multiplication.^[5][6]

This section has a simplified version of the algorithm, showing how to compute the product${\displaystyle ab}$ of two natural numbers${\displaystyle a,b}$, modulo a number of the form${\displaystyle 2^n+1}$, where${\displaystyle n=2^{k}M}$ is some fixed number. The integers${\displaystyle a,b}$ are to be divided into${\displaystyle D=2^k}$ blocks of${\displaystyle M}$ bits, so in practical implementations, it is important to strike the right balance between the parameters${\displaystyle M,k}$. In any case, this algorithm will provide a way to multiply two positive integers, provided${\displaystyle n}$ is chosen so that.

Let${\displaystyle n=DM}$ be the number of bits in the signals${\displaystyle a}$ and${\displaystyle b}$, where${\displaystyle D=2^k}$ is a power of two. Divide the signals${\displaystyle a}$ and${\displaystyle b}$ into${\displaystyle D}$ blocks of${\displaystyle M}$ bits each, storing the resulting blocks as arrays${\displaystyle A,B}$ (whose entries we shall consider for simplicity as arbitrary precision integers).

We now select a modulus for the Fourier transform, as follows. Let${\displaystyle M^{\prime }}$ be such that${\displaystyle D{M}^{\prime }\ge 2M+k}$. Also put${\displaystyle n^{\prime }=D{M}^{\prime }}$, and regard the elements of the arrays${\displaystyle A,B}$ as (arbitrary precision) integers modulo${\displaystyle 2^{n^{\prime }}+1}$. Observe that since${\displaystyle 2^{n^{\prime }}+1\ge 2^{2M+k}+1=D{2}^{2M}+1}$, the modulus is large enough to accommodate any carries that can result from multiplying${\displaystyle a}$ and${\displaystyle b}$. Thus, the product${\displaystyle ab}$ (modulo${\displaystyle 2^n+1}$) can be calculated by evaluating the convolution of${\displaystyle A,B}$. Also, with${\displaystyle g=2^{2M^{\prime }}}$, we have${\displaystyle g^{D/2}\equiv -1(\mathrm{mod}{2}^{n^{\prime }}+1)}$, and so${\displaystyle g}$ is a primitive${\displaystyle D}$th root of unity modulo${\displaystyle 2^{n^{\prime }}+1}$.

We now take the discrete Fourier transform of the arrays${\displaystyle A,B}$ in the ring${\displaystyle \mathbb{Z}/(2^{n^{\prime }}+1)\mathbb{Z}}$, using the root of unity${\displaystyle g}$ for the Fourier basis, giving the transformed arrays${\displaystyle \hat{A},\hat{B}}$. Because${\displaystyle D=2^k}$ is a power of two, this can be achieved in logarithmic time using afast Fourier transform.

Let${\displaystyle {\hat{C}}_i={\hat{A}}_i{\hat{B}}_i}$ (pointwise product), and compute the inverse transform${\displaystyle C}$ of the array${\displaystyle \hat{C}}$, again using the root of unity${\displaystyle g}$. The array${\displaystyle C}$ is now the convolution of the arrays${\displaystyle A,B}$. Finally, the product${\displaystyle ab(\mathrm{mod}{2}^n+1)}$ is given by evaluating$${\displaystyle ab\equiv \sum _j{C}_j{2}^{Mj}\mathrm{mod}{2}^n+1.}$$

This basic algorithm can be improved in several ways. Firstly, it is not necessary to store the digits of${\displaystyle a,b}$ to arbitrary precision, but rather only up to${\displaystyle n^{\prime }+1}$ bits, which gives a more efficient machine representation of the arrays${\displaystyle A,B}$. Secondly, it is clear that the multiplications in the forward transforms are simple bit shifts. With some care, it is also possible to compute the inverse transform using only shifts. Taking care, it is thus possible to eliminate any true multiplications from the algorithm except for where the pointwise product${\displaystyle {\hat{C}}_i={\hat{A}}_i{\hat{B}}_i}$ is evaluated. It is therefore advantageous to select the parameters${\displaystyle D,M}$ so that this pointwise product can be performed efficiently, either because it is a single machine word or using some optimized algorithm for multiplying integers of a (ideally small) number of words. Selecting the parameters${\displaystyle D,M}$ is thus an important area for further optimization of the method.

Every number in base B, can be written as a polynomial:

${\displaystyle X=\sum _{i=0}^N{x}_i{B}^i}$

Furthermore, multiplication of two numbers could be thought of as a product of two polynomials:

${\displaystyle XY=\left(\sum _{i=0}^N{x}_i{B}^i\right)\left(\sum _{j=0}^N{y}_i{B}^j\right)}$

Because, for${\displaystyle B^k}$:${\displaystyle c_k=\sum _{(i,j):i+j=k}{a}_i{b}_j=\sum _{i=0}^k{a}_i{b}_{k-i}}$,we have a convolution.

By using FFT (fast Fourier transform), used in the original version rather than NTT (Number-theoretic transform),^[7] with convolution rule; we get

${\displaystyle \hat{f}(a\ast b)=\hat{f}\left(\sum _{i=0}^k{a}_i{b}_{k-i}\right)=\hat{f}(a)\bullet \hat{f}(b).}$

That is;${\displaystyle C_k=a_k\bullet b_k}$, where${\displaystyle C_k}$is the corresponding coefficient in Fourier space. This can also be written as:${\displaystyle \text{fft}(a\ast b)=\text{fft}(a)\bullet \text{fft}(b)}$.

We have the same coefficients due to linearity under the Fourier transform, and because these polynomials only consist of one unique term per coefficient:

${\displaystyle \hat{f}(x^n)={\left(\frac{i}{2\pi }\right)}^n{\delta }^{(n)}}$ and

${\displaystyle \hat{f}(aX(\xi )+bY(\xi ))=a\hat{X}(\xi )+b\hat{Y}(\xi )}$

Convolution rule:${\displaystyle \hat{f}(X\ast Y)=\hat{f}(X)\bullet \hat{f}(Y)}$

We have reduced our convolution problemto product problem, through FFT.

By finding the FFT of thepolynomial interpolation of each${\displaystyle C_k}$, one can determine the desired coefficients.

This algorithm uses thedivide-and-conquer method to divide the problem into subproblems.

${\displaystyle c_k=\sum _{(i,j):i+j\equiv k(\mathrm{mod}N(n))}{a}_i{b}_j}$, where${\displaystyle N(n)=2^n+1}$.

By letting:

${\displaystyle a_i^{\prime }={\theta }^i{a}_i}$ and${\displaystyle b_j^{\prime }={\theta }^j{b}_j,}$

where${\displaystyle {\theta }^N=-1}$ is the n^th root, one sees that:^[8]

This mean, one can use weight${\displaystyle {\theta }^i}$, and then multiply with${\displaystyle {\theta }^{-k}}$ after.

Instead of using weight, as${\displaystyle {\theta }^N=-1}$, in first step of recursion (when${\displaystyle n=N}$), one can calculate:

${\displaystyle C_k=\sum _{(i,j):i+j\equiv k(\mathrm{mod}N(N))}=\sum _{(i,j):i+j=k}{a}_i{b}_j-\sum _{(i,j):i+j=k+n}{a}_i{b}_j}$

In a normal FFT which operates over complex numbers, one would use:

$${\displaystyle \exp \left(\frac{2k\pi i}{n}\right)=\cos \frac{2k\pi }{n}+i\sin \frac{2k\pi }{n},k=0,1,\dots ,n-1.}$$

However, FFT can also be used as a NTT (number theoretic transformation) in Schönhage–Strassen. This means that we have to use θ to generate numbers in a finite field (for example${\displaystyle \mathrm{G}\mathrm{F}(2^n+1)}$).

A root of unity under a finite fieldGF(r), is an element a such that${\displaystyle {\theta }^{r-1}\equiv 1}$ or${\displaystyle {\theta }^r\equiv \theta }$. For exampleGF(p), wherep is aprime number, gives${\displaystyle \{1,2,\dots ,p-1\}}$.

Notice that${\displaystyle 2^n\equiv -1}$ in${\displaystyle \mathrm{GF}(2^n+1)}$ and${\displaystyle \sqrt{2}\equiv -1}$ in${\displaystyle \mathrm{GF}(2^{n+2}+1)}$. For these candidates,${\displaystyle {\theta }^N\equiv -1}$ under its finite field, and therefore act the way we want .

Same FFT algorithms can still be used, though, as long asθ is aroot of unity of a finite field.

To find FFT/NTT transform, we do the following:

First product gives contribution to${\displaystyle c_k}$, for eachk. Second gives contribution to${\displaystyle c_k}$, due to${\displaystyle (i+j)}$ mod${\displaystyle N(n)}$.

To do the inverse:

${\displaystyle C_k=2^{-m}\widehat{f^{-1}}({\theta }^{-k}{C}_{k+k}^{\prime })}$ or${\displaystyle C_k=\widehat{f^{-1}}({\theta }^{-k}{C}_{k+k}^{\prime })}$

depending whether data needs to be normalized.

One multiplies by${\displaystyle 2^{-m}}$ to normalize FFT data into a specific range, where${\displaystyle \frac{1}{n}\equiv 2^{-m}modN(n)}$, wherem is found using themodular multiplicative inverse.

In Schönhage–Strassen algorithm,${\displaystyle N=2^M+1}$. This should be thought of as a binary tree, where one have values in${\displaystyle 0\le \text{index}\le 2^M=2^{i+j}}$. By letting${\displaystyle K\in [0,M]}$, for eachK one can find all${\displaystyle i+j=K}$, and group all${\displaystyle (i,j)}$ pairs into M different groups. Using${\displaystyle i+j=k}$ to group${\displaystyle (i,j)}$ pairs through convolution is a classical problem in algorithms.^[9]

Having this in mind,${\displaystyle N=2^M+1}$ help us to group${\displaystyle (i,j)}$ into${\displaystyle \frac{M}{2^k}}$ groups for each group of subtasks in depthk in a tree with${\displaystyle N=2^{\frac{M}{2^k}}+1}$

Notice that${\displaystyle N=2^M+1=2^{2^L}+1}$, for some L. This makes N aFermat number. When doing mod${\displaystyle N=2^M+1=2^{2^L}+1}$, we have a Fermat ring.

Because some Fermat numbers are Fermat primes, one can in some cases avoid calculations.

There are otherN that could have been used, of course, with same prime number advantages. By letting${\displaystyle N=2^k-1}$, one have the maximal number in a binary number with${\displaystyle k+1}$ bits.${\displaystyle N=2^k-1}$ is a Mersenne number, that in some cases is a Mersenne prime. It is a natural candidate against Fermat number${\displaystyle N=2^{2^L}+1}$

Doing several mod calculations against differentN, can be helpful when it comes to solving integer product. By using theChinese remainder theorem, after splittingM into smaller different types ofN, one can find the answer of multiplicationxy^[10]

Fermat numbers and Mersenne numbers are just two types of numbers, in something called generalized Fermat Mersenne number (GSM); with formula:^[11]

${\displaystyle G_{q,p,n}=\sum _{i=1}^p{q}^{(p-i)n}=\frac{q^{pn}-1}{q^n-1}}$

${\displaystyle M_{p,n}=G_{2,p,n}}$

In this formula,${\displaystyle M_{2,2^k}}$ is a Fermat number, and${\displaystyle M_{p,1}}$ is a Mersenne number.

This formula can be used to generate sets of equations, that can be used in CRT (Chinese remainder theorem):^[12]

${\displaystyle g^{\frac{(M_{p,n}-1)}{2}}\equiv -1(\mathrm{mod}{M}_{p,n})}$, whereg is a number such that there exists anx where${\displaystyle x^2\equiv g(\mathrm{mod}{M}_{p,n})}$, assuming${\displaystyle N=2^n}$

Furthermore;${\displaystyle g^{2^{(p-1)n}-1}\equiv a^{2^n-1}(\mathrm{mod}{M}_{p,n})}$, wherea is an element that generates elements in${\displaystyle \{1,2,4,..{.2}^{n-1},2^n\}}$ in a cyclic manner.

If${\displaystyle N=2^t}$, where${\displaystyle 1\le t\le n}$, then${\displaystyle g_t=a^{(2^n-1)2^{n-t}}}$.

The following formula is helpful, finding a properK (number of groups to divideN bits into) given bit sizeN by calculating efficiency :^[13]

${\displaystyle E=\frac{\frac{2N}{K}+k}{n}}$N is bit size (the one used in${\displaystyle 2^N+1}$) at outermost level.K gives${\displaystyle \frac{N}{K}}$ groups of bits, where${\displaystyle K=2^k}$.

n is found throughN, K andk by finding the smallestx, such that${\displaystyle 2N/K+k\le n=K{2}^x}$

If one assume efficiency above 50%,${\displaystyle \frac{n}{2}\le \frac{2N}{K},K\le n}$ andk is very small compared to rest of formula; one get

${\displaystyle K\le 2\sqrt{N}}$

This means: When something is very effective;K is bound above by${\displaystyle 2\sqrt{N}}$ or asymptotically bound above by${\displaystyle \sqrt{N}}$

Following algorithm, the standard Modular Schönhage-Strassen Multiplication algorithm (with some optimizations), is found in overview through^[14]

• T3MUL = Toom–Cook multiplication
• SMUL = Schönhage–Strassen multiplication
• Evaluate = FFT/IFFT

For implementation details, one can read the bookPrime Numbers: A Computational Perspective.^[15] This variant differs somewhat from Schönhage's original method in that it exploits thediscrete weighted transform to performnegacyclic convolutions more efficiently. Another source for detailed information isKnuth'sThe Art of Computer Programming.^[16]

This section explains a number of important practical optimizations, when implementing Schönhage–Strassen.

Below a certain cutoff point, it's more efficient to use other multiplication algorithms, such asToom–Cook multiplication.^[17]

The idea is to use${\displaystyle \sqrt{2}}$ as aroot of unity of order${\displaystyle 2^{n+2}}$ in finite field${\displaystyle \mathrm{G}\mathrm{F}(2^{n+2}+1)}$ (it is a solution to equation${\displaystyle {\theta }^{2^{n+2}}\equiv 1(\mathrm{mod}{2}^{n+2}+1)}$), when weighting values in NTT (number theoretic transformation) approach. It has been shown to save 10% in integer multiplication time.^[18]

By letting${\displaystyle m=N+h}$, one can compute${\displaystyle uvmod{2}^N+1}$ and${\displaystyle (umod{2}^h)(v{mod2}^h)}$ in combination with CRT (Chinese Remainder Theorem) to find exact values of multiplicationuv.^[19]

• Computer arithmetic algorithms
• Multiplication

• CS1 German-language sources (de)
• Webarchive template wayback links
• CS1: long volume value
• Articles with short description
• Short description matches Wikidata
• Articles needing additional references from October 2024
• All articles needing additional references