Schönhage–Strassen algorithm

Details

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), 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 problem to product problem, through FFT.

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

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

Convolution under mod N

${\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:

${\displaystyle \begin{array}{cc}{C}_k& =\sum _{(i,j):i+j\equiv k(\mathrm{mod}N(n))}{a}_i{b}_j={\theta }^{-k}\sum _{(i,j):i+j\equiv k(\mathrm{mod}N(n))}{a}_i^{\prime }{b}_j^{\prime }\\ & ={\theta }^{-k}\left(\sum _{(i,j):i+j=k}{a}_i^{\prime }{b}_j^{\prime }+\sum _{(i,j):i+j=k+n}{a}_i^{\prime }{b}_j^{\prime }\right)\\ & ={\theta }^{-k}\left(\sum _{(i,j):i+j=k}{a}_i{b}_j{\theta }^k+\sum _{(i,j):i+j=k+n}{a}_i{b}_j{\theta }^{n+k}\right)\\ & =\sum _{(i,j):i+j=k}{a}_i{b}_j+{\theta }^n\sum _{(i,j):i+j=k+n}{a}_i{b}_j.\end{array}}$

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.}$$

${\displaystyle \begin{array}{cc}{C}_k& ={\theta }^{-k}\left(\sum _{(i,j):i+j=k}{a}_i{b}_j{\theta }^k+\sum _{(i,j):i+j=k+n}{a}_i{b}_j{\theta }^{n+k}\right)\\ & =e^{-i2\pi k/n}\left(\sum _{(i,j):i+j=k}{a}_i{b}_j{e}^{i2\pi k/n}+\sum _{(i,j):i+j=k+n}{a}_i{b}_j{e}^{i2\pi (n+k)/n}\right)\end{array}}$

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 field GF(r), is an element a such that ${\displaystyle {\theta }^{r-1}\equiv 1}$ or ${\displaystyle {\theta }^r\equiv \theta }$. For example GF(p), where p is a prime 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 a root of unity of a finite field.

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

${\displaystyle \begin{array}{cc}{C}_k^{\prime }& =\hat{f}(k)=\hat{f}\left({\theta }^{-k}\left(\sum _{(i,j):i+j=k}{a}_i{b}_j{\theta }^k+\sum _{(i,j):i+j=k+n}{a}_i{b}_j{\theta }^{n+k}\right)\right)\\ C_{k+k}^{\prime }& =\hat{f}(k+k)=\hat{f}\left(\sum _{(i,j):i+j=2k}{a}_i{b}_j{\theta }^k+\sum _{(i,j):i+j=n+2k}{a}_i{b}_j{\theta }^{n+k}\right)\\ & =\hat{f}\left(\sum _{(i,j):i+j=2k}{a}_i{b}_j{\theta }^k+\sum _{(i,j):i+j=2k+n}{a}_i{b}_j{\theta }^{n+k}\right)\\ & =\hat{f}\left(A_{k\leftarrow k}\right)\bullet \hat{f}(B_{k\leftarrow k})+\hat{f}(A_{k\leftarrow k+n})\bullet \hat{f}(B_{k\leftarrow k+n})\end{array}}$

First product gives contribution to ${\displaystyle c_k}$, for each k. 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)}$, where m is found using the modular multiplicative inverse.

Implementation details

Why N = 2M + 1 mod N

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 each K 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.

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 depth k 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 a Fermat 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 other N 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}$

In search of another N

Doing several mod calculations against different N, can be helpful when it comes to solving integer product. By using the Chinese remainder theorem, after splitting M into smaller different types of N, one can find the answer of multiplication xy

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

${\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):

${\displaystyle g^{\frac{(M_{p,n}-1)}{2}}\equiv -1(\mathrm{mod}{M}_{p,n})}$, where g is a number such that there exists an x 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})}$, where a 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}}}$.

How to choose K for a specific N

The following formula is helpful, finding a proper K (number of groups to divide N bits into) given bit size N by calculating efficiency :

${\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 through N, K and k by finding the smallest x, 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}$ and k 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}}$

Pseudocode

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

1. Split both input numbers a and b into n coefficients of s bits each. Use at least ⁠${\displaystyle K+1}$⁠ bits to store them, to allow encoding of the value ⁠${\displaystyle 2^K.}$⁠
2. Weight both coefficient vectors according to (2.24) with powers of θ by performing cyclic shifts on them.
3. Shuffle the coefficients ⁠${\displaystyle a_i}$⁠ and ⁠${\displaystyle b_j}$⁠ .
4. Evaluate ⁠${\displaystyle a_i}$⁠ and ⁠${\displaystyle b_j}$⁠ . Multiplications by powers of ω are cyclic shifts.
5. Do n pointwise multiplications ⁠${\displaystyle c_k:=a_k{b}_k}$⁠ in ⁠${\displaystyle Z/(2^K+1)Z}$⁠. If SMUL is used recursively, provide K as parameter. Otherwise, use some other multiplication function like T3MUL and reduce modulo ⁠${\displaystyle 2^K+1}$⁠ afterwards.
6. Shuffle the product coefficients ⁠${\displaystyle c_k}$⁠.
7. Evaluate the product coefficients ⁠${\displaystyle c_k}$⁠.
8. Apply the counterweights to the ⁠${\displaystyle c_k}$⁠ according to (2.25). Since ⁠${\displaystyle {\theta }^{2n}\equiv 1}$⁠ it follows that ⁠${\displaystyle {\theta }^{-k}\equiv {\theta }^{n-k}}$⁠
9. Normalize the ⁠${\displaystyle c_k}$⁠ with ⁠${\displaystyle 1/n\equiv 2^{-m}}$⁠ (again a cyclic shift).
10. Add up the ⁠${\displaystyle c_k}$⁠ and propagate the carries. Make sure to properly handle negative coefficients.
11. Do a reduction modulo ⁠${\displaystyle 2^N+1}$⁠.

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

Further study

For implementation details, one can read the book Prime Numbers: A Computational Perspective. 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.

Optimizations

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

Use of other multiplications algorithm, inside algorithm

Below a certain cutoff point, it's more efficient to use other multiplication algorithms, such as Toom–Cook multiplication.

Square root of 2 trick

The idea is to use ${\displaystyle \sqrt{2}}$ as a root 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.

Granlund's trick

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 multiplication uv.