Convolution under mod N

${\displaystyle c_{k}=\sum _{(i,j):i+j\equiv k{\pmod {N(n)}}}a_{i}b_{j}}$, where ${\displaystyle N(n)=2^{n}+1}$ and ${\displaystyle N(N)=2^{N}+1}$ in Schönhage–Strassen algorithm.

By letting:

${\displaystyle a_{i}'=\theta ^{i}a_{i}}$ and ${\displaystyle b_{j}'=\theta ^{j}b_{j},}$

where ${\displaystyle \theta ^{N}=-1}$ is the n-th root

One sees that:^[8]

${\displaystyle {\begin{aligned}C_{k}&=\sum _{(i,j):i+j=k\equiv {\pmod {N(n)}}}a_{i}b_{j}=\theta ^{-k}\sum _{(i,j):i+j\equiv k{\pmod {N(n)}}}a_{i}'b_{j}'\\[6pt]&=\theta ^{-k}\left(\sum _{(i,j):i+j=k}a_{i}'b_{j}'+\sum _{(i,j):i+j=k+n}a_{i}'b_{j}'\right)\\[6pt]&=\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)\\[6pt]&=\sum _{(i,j):i+j=k}a_{i}b_{j}+\theta ^{n}\sum _{(i,j):i+j=k+n}a_{i}b_{j}.\end{aligned}}}$

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

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

${\displaystyle C_{k}=\sum _{(i,j):i+j\equiv k{\pmod {N(N)}}}=\sum _{(i,j):i+j=k}a_{i}b_{j}-\sum _{(i,j):i+j=k+n}a_{i}b_{j}}$

In normal FFT, that 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}},\qquad k=0,1,\dots ,n-1.}$$

${\displaystyle {\begin{aligned}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)\\[6pt]&=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{aligned}}}$

However, FFT can also be used as a NTT (number theoretic transformation) in Schönhage–Strassen. This means that we have to use θ that generate numbers in a finite field (for example ${\displaystyle \mathrm {GF} (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, gives ${\displaystyle \{1,2,\ldots ,p-1\}}$.

Notice that ${\displaystyle 2^{n}\equiv -1}$ in ${\displaystyle \operatorname {GF} (2^{n}+1)}$ and ${\displaystyle {\sqrt {2}}\equiv -1}$ in ${\displaystyle \operatorname {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 root of unity of a finite field.

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

${\displaystyle {\begin{aligned}C_{k}'&={\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)\\[6pt]C_{k+k}'&={\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)\\[6pt]&={\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)\\[6pt]&={\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{aligned}}}$

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}{\hat {f^{-1}}}(\theta ^{-k}C_{k+k}')}$ or ${\displaystyle C_{k}={\hat {f^{-1}}}(\theta ^{-k}C_{k+k}')}$

depending on whether fft one use normalize data or not.

One multiply by ${\displaystyle 2^{-m}}$, to normailize fft data to a specific range, where ${\displaystyle {\frac {1}{n}}\equiv 2^{-m}{\bmod {N}}(n)}$, where m is found using modular multiplicative inverse.

Why N = 2^M + 1 in mod N

In Schönhage–Strassen algorithm, ${\displaystyle N=2^{M}+1}$. One should think of this as a binary tree, where one have values in ${\displaystyle 0\leq {\text{index}}\leq 2^{M}=2^{i+j}}$. By letting ${\displaystyle K\in [0,M]}$, one can for each K find all ${\displaystyle i+j=K}$: One can 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] For example: Let k be total income and i be mans income and j womans income; by using convolution, one can group ${\displaystyle (i,j)}$ into K groups based on desired total income.

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 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 is Fermat number. When doing mod ${\displaystyle N=2^{M}+1=2^{2^{L}}+1}$, we have something called 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}$

How to choose K for a specific N

Following formula helps one out, finding a proper K (number of groups to divide N bits into) given bit size N 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 through N, K and k by finding the smallest x, such that ${\displaystyle 2N/K+k\leq n=K2^{x}}$

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

${\displaystyle K\leq 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 alogithm, 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

Further study

For implemantion details, one can read the book Prime Numbers: A Computational Perspective.^[15] 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.^[16]

Use of other multiplications algorithm, inside algorithm

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

Square root of 2 trick

The idea is to use ${\displaystyle {\sqrt {2}}}$ as root of unity of order ${\displaystyle 2^{n+2}}$ in finite field ${\displaystyle \mathrm {GF} (2^{n+2}+1)}$ ( it is a solution to equation ${\displaystyle \theta ^{2^{n+2}}\equiv 1{\pmod {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]

Granlund's trick

By letting ${\displaystyle m=N+h}$; one can compute ${\displaystyle uv{\bmod {2^{N}+1}}}$ and ${\displaystyle (u{\bmod {2^{h}}})(v{\bmod {2}}^{h})}$. In combination with CRT (Chinese Remainder Theorem), to find exact values of multiplication uv^[19]