Number Theoretic Transform

Number Theoretic Transform(NTT) is an efficient algorithm for polynomial multiplication, akin to the Fast Fourier Transform (FFT). It operates in finite fields, making it particularly valuable in cryptography and coding theory, where it enables rapid computations on large numbers and polynomials. For the purposes of this document, polynomials are assumed to be defined over finite fields by default.

Polynomial Multiplication

Let's consider two polynomials: $A(x)=x^2+3x+2$        $B(x)=2x^2+1$
The product $C(x)=A(x)\cdot B(x)$ can be expressed as:
$C(x)=(x^2+3x+2)(2x^2+1)$
By expanding, we find: $C(x)=x^2(2x^2+1)+3x(2x^2+1)+2(2x^2+1)$
$C(x)=2x^4+x^2+6x^3+3x+4x^2+2$
$C(x)=2x^4+6x^3+5x^2+3x+2$

In a programming context, we can represent polynomials as follows:
$A(x)=2+3x+x^2\to A=[2,3,1]$
$B(x)=1+2x^2\to B=[1,0,2]$
$C(x)=2+3x+5x^2+6x^3+2x^4\to C=[2,3,5,6,2]$
$C[k]$ = coefficient of $k^{th}$ degree term of polynomial $C(x)$

$A(x)=a_0+a_{1}x+a_2{x}^2+...+a_d{x}^d$
$B(x)=b_0+b_{1}x+b_2{x}^2+...+b_d{x}^d$
$C(x)=c_0+c_{1}x+c_2{x}^2+...+c_{2d}{x}^{2d}$
Using the distributive property, the runtime for multiplying two polynomials of degree $d$ is $O(d^2)$.However, is there a more efficient method?

Polynomial Representation

Coefficient Representation

$P(x)=p_0+p_{1}x\to [p_0,p_1]$
To illustrate this representation visually,we can use an interactive graph where the coefficients can be adjusted:

Value Representation

$P(x)=p_0+p_{1}x\to [(x_0,P(x_0)),(x_1,P(x_1))]$ A polynomial can also be uniquely defined by its values at some specific points. For example, two points define a unique line:

• For points$(-2,0)$ and $(2,2)$, we have $P(x)=1+0.5s$.
• For points$(3,0)$ and $(-1,4)$, we have $P(x)=3-x$.

Three points define a quadratic polynomial:

• $\{(-3,1),(-1,-1),(1,3)\}\to P(x)=\frac{3}{4}{x}^2+2x+\frac{1}{4}$
• $\{(-1,0),(0,1),(1,0),(2,1)\}\to P(x)=\frac{2}{3}{x}^3-x^2-\frac{2}{3}x+1$

More generally,$(d+1)$ points uniquely define a polynomial of degree $d$: $$\{(x_0,P(x_0)),(x_1,P(x_1)),...,(x_d,P(x_d))\}$$ $$P(x)=p_0+p_{1}x+p_2{x}^2+...+p_d{x}^d$$ $$\begin{gathered} P(x_0)=p_0+p_1{x}_0+p_2{x}_0^2+...+p_d{x}_0^d \\ P(x_1)=p_0+p_1{x}_1+p_2{x}_1^2+...+p_d{x}_1^d \\ ⋮ \\ P(x_d)=p_0+p_1{x}_d+p_2{x}_d^2+...+p_d{x}_d^d \end{gathered}$$ That is: $$\begin{array}{c}\left[\begin{array}{c}P(x_0)\\ P(x_1)\\ ⋮\\ P(x_d)\end{array}\right]=\underset{M}{\underbrace{\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^d\\ 1& x_1& x_1^2& \cdots & x_1^d\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_d& x_d^2& \cdots & x_d^d\end{array}\right]}}\left[\begin{array}{c}{p}_0\\ p_1\\ ⋮\\ p_d\end{array}\right]\end{array}$$ $M$ is invertible for unique $x_0,x_1,...,x_d$ $\Rightarrow$ Unique $p_0,p_1,...,p_d$ exists $\Rightarrow$ Unique polynomial $P(x)$ exists.

Punchline: Two Unique Trpresentations for Polynomials $$P(x)=p_0+p_{1}x+p_2{x}^2+...+p_d{x}^d$$ $$1.\underset{\text{Coefficient Representation}}{\underbrace{[p_0,p_1,...,p_d]}}$$ $$2.\underset{\text{Value Representation}}{\underbrace{(x_0,P(x_0),(x_1,P(x_1)),...,(x_d,P(x_d)))}}$$

Value Representation Advantages

$\underset{[(-2,1),(-1,0),(0,1),(1,4),(2,9)]}{A(x)=x^2+2x+1}\underset{[(-2,9),(-1,4),(0,1),(1,0),(2,1)]}{B(x)=x^2-2x+1}\underset{Degree4\to 5points}{\underbrace{C(x)=A(x)\cdot B(x)}}$
and $\begin{array}{r}[(-2,1),(-1,0),(0,1),(1,4),(2,9)]\\ \times [(-2,9),(-1,4),(0,1),(1,0),(2,1)]\\ [(-2,9),(-1,0),(0,1),(1,0),(2,9)]\end{array}$
$\Rightarrow$ $\underset{[(-2,9),(-1,0),(0,1),(1,0),(2,9)]}{C(x)=A(x)\cdot B(x)}$ $\Rightarrow$ $\underset{[(-2,9),(-1,0),(0,1),(1,0),(2,9)]}{C(x)=x^4-2x^2+1}$

In value representation, multiplication can be performed in linear time $O(d)$ rather than quadratic $O(d^2)$. This is especially advantageous when working with large polynomials.

Polynomial Multiplication Flowchart

Polynomial Evaluation

Evaluation from Coefficients to Values

To evaluate a polynomial $P(x)=p_0+p_{1}x+p_2{x}^2+...+p_d{x}^d$ at $n\ge d+1$ points,we can compute as follows: $$\begin{array}{c}(1,P(1))\\ (2,P(2))\\ ⋮\\ (n,P(n))\end{array}\Rightarrow \underset{O(nd)⟺O(d^2)}{\underbrace{\begin{array}{c}(1,p_0+p_1\cdot 1+p_2\cdot 1^2+\cdots +p_d\cdot 1^d)\\ (2,p_0+p_1\cdot 2+p_2\cdot 2^2+\cdots +p_d\cdot 2^d)\\ ⋮\\ (n,p_0+p_1\cdot n+p_2\cdot n^2+\cdots +p_d\cdot n^d)\end{array}}}$$ Howerer, this operation can have a complexity of $O(nd)$, which could lead to $O(d^2)$.

More Efficient Strategies

For example, when evaluatiing $P(x)=x^2$ at 8 points, we can select symmetric points such as
$$\begin{gathered} (1,1)(-1,1) \\ (2,4)(-2,4) \\ (3,9)(-3,9) \\ (4,16)(-4,16) \end{gathered}$$.
We have $P(-x)=P(x)$.

For $P(x)=x^3$ we can select:
$$\begin{gathered} (1,1)(-1,-1) \\ (2,8)(-2,-8) \\ (3,27)(-3,-27) \\ (4,64)(-4,-64) \end{gathered}$$
We have $P(-x)=-P(x)$.

So we only need 4 points!

• $$P(x)=3x^5+2x^4+x^3+7x^2+5x+1$$ $$Evaluateatnpoints\pm x_1,\pm x_2,\dots ,\pm x_{\frac{n}{2}}$$ $$P(x)=(2x^4+7x^2+1)+(3x^5+x^3+5x)=\underset{P_e(x^2)}{\underbrace{(2x^4+7x^2+1)}}+x\underset{P_o(x^2)}{\underbrace{(3x^4+x^2+5)}}$$ $$P(x)=P_e(x^2)+x{P}_o(x^2)$$ $$\begin{array}{l}P(x_i)=P_e(x_i^2)+x_i{P}_o(x_i^2)\\ P(-x_i)=P_e(x_i^2)-x_i{P}_o(x_i^2)\end{array}\}\text{Lot of overlap!}$$ $$P_e(x^2)=2x^2+7x+1P_o(x^2)=3x^2+x+5$$ $$P_e(x^2)\text{ and }{P}_o(x^2)\text{ have degree 2!}$$

• $$P(x)=p_0+p_{1}x+p_2{x}^2+\cdots +p_{n-1}{x}^{n-1}$$ $$\text{Evaluate at }n\text{ points }\pm x_1,\pm x_2,\dots ,\pm x_{\frac{n}{2}}$$ $$P(x)=P_e(x^2)+x{P}_o(x^2)$$ $$\begin{array}{l}P(x_i)=P_e(x_i^2)+x_i{P}_o(x_i^2)\\ P(-x_i)=P_e(x_i^2)-x_i{P}_o(x_i^2)\end{array}\}\text{Lot of overlap!}$$ $$P_e(x^2)\text{ and }{P}_o(x^2)\text{ have degree }\frac{n}{2}!$$

$$\underset{\text{Same process on simpler problem}}{\underbrace{\text{Evaluate }{P}_e(x^2)\text{ and }{P}_o(x^2)\text{ each at }{x}_1^2,x_2^2,\dots ,x_{\frac{n}{2}}^2(\frac{n}{2}\text{ points})}}$$

But that is one major problem. $$\text{Points }[\pm x_1,\pm x_2,\dots ,\pm x_{fracn2}]are\pm paired.$$ $$\text{Points }[x_1^2,x_2^2,\dots ,x_{\frac{n}{2}}^2]\text{ are not }\pm \text{ paired.}$$ $$\text{Recursion breaks!}$$ $$\text{Is it possible to make }[x_1^2,x_2^2,\dots ,x_{\frac{n}{2}}^2\pm \text{ paired?}$$ We can choose the roots of unity in the finite field.

Which Evaluation Points to Use?

To evaluation $P(x)=p_0+p_{1}x+p_2{x}^2+\cdots +p_d{x}^d$, we need at least $n\ge (d+1)$ points, ideally $n=2^k$ for $k\in \mathbb{Z}$. The points chosen should be the $n^{th}$ roots of unity:

$$\text{Evaluate }P(x)\text{ at }[1,{\omega }^1,{\omega }^2,\dots ,{\omega }^{n-1}]$$ whers $\omega$ is a root of $z^n=1$.

NTT Implementation

def FFT(P):
    # P- [p0,p1,...,p_{n-1}] coeff representation
    n=len(P) # n is a power of 2
    if n==1:
        return P
    w = get_root_of_unity(n)
    P_e,P_o = [P[0],P[2],...,P[n-2]],[P[1],P[3],...,P[n-1]]
    y_e,y_o = FFT(P_e), FFT(P_o)
    y = [0]*n
    for j in range(n/2):
        y[j] = y_e[j]+w^j*y_o[j]
        y[j+n/2] = y_e[j]-w^j*y_o[j]
    return y

Interpolation and Inverse NTT

Alternative Perspective on Evaluation/NTT

$$P(x)=p_0+p_{1}x+p_2{x}^2+\cdots +p_{n-1}{x}^{n-1}$$ $$\begin{array}{c}\left[\begin{array}{c}P(x_0)\\ P(x_1)\\ ⋮\\ P(x_{n-1})\end{array}\right]=\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^{n-1}\\ 1& x_1& x_1^2& \cdots & x_1^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n-1}& x_{n-1}^2& \cdots & x_{n-1}^{n-1}\end{array}\right]\left[\begin{array}{c}{p}_0\\ p_1\\ ⋮\\ p_{n-1}\end{array}\right]\end{array}$$
$$x_k={\omega }^k\text{ where }\omega \text{ is a root of unity in the finite field}$$ That is $$\begin{array}{c}\left[\begin{array}{c}P(x_0)\\ P(x_1)\\ ⋮\\ P(x_{n-1})\end{array}\right]=\underset{\text{Discrete Fourier Transform (DFT) matrix}}{\underbrace{\left[\begin{array}{ccccc}1& 1& 1^2& \cdots & 1^{n-1}\\ 1& \omega & {\omega }^2& \cdots & {\omega }^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }^{n-1}& {\omega }^{2(n-1)}& \cdots & {\omega }^{(n-1)(n-1)}\end{array}\right]}}\left[\begin{array}{c}{p}_0\\ p_1\\ ⋮\\ p_{n-1}\end{array}\right]\end{array}$$

Interpolation involves inversing the DFT matrix

$$P(x)=p_0+p_{1}x+p_2{x}^2+\cdots +p_{n-1}{x}^{n-1}$$ $$\begin{array}{c}\begin{array}{c}\left[\begin{array}{c}{p}_0\\ p_1\\ ⋮\\ p_{n-1}\end{array}\right]={\left[\begin{array}{ccccc}1& 1& 1^2& \cdots & 1^{n-1}\\ 1& \omega & {\omega }^2& \cdots & {\omega }^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }^{n-1}& {\omega }^{2(n-1)}& \cdots & {\omega }^{(n-1)(n-1)}\end{array}\right]}^{-1}\left[\begin{array}{c}P(x_0)\\ P(x_1)\\ ⋮\\ P(x_{n-1})\end{array}\right]\end{array}\\ \Downarrow \\ \begin{array}{c}\left[\begin{array}{c}{p}_0\\ p_1\\ ⋮\\ p_{n-1}\end{array}\right]=\frac{1}{n}\left[\begin{array}{ccccc}1& 1& 1^2& \cdots & 1^{n-1}\\ 1& {\omega }^{-1}& {\omega }^{-2}& \cdots & {\omega }^{-(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }^{-(n-1)}& {\omega }^{-2(n-1)}& \cdots & {\omega }^{-(n-1)(n-1)}\end{array}\right]\left[\begin{array}{c}P(x_0)\\ P(x_1)\\ ⋮\\ P(x_{n-1})\end{array}\right]\end{array}\end{array}$$
The inverse matrix and original matrix look quite similar! Every $\omega$ in original matrix is now $\frac{1}{n}{\omega }^{-1}$

$$INTT(<values>)⟺NTT(<values>)\text{ wit }{\omega }^{\prime }=\frac{1}{n}{\omega }^{-1}$$

So the implementation of inverse NTT as follows:

def IFFT(P):
    # P- [p0,p1,...,p_{n-1}] coeff representation
    n=len(P) # n is a power of 2
    if n==1:
        return P
    w = 1/n*(get_root_of_unity(n)^{-1})
    P_e,P_o = [P[0],P[2],...,P[n-2]],[P[1],P[3],...,P[n-1]]
    y_e,y_o = FFT(P_e), FFT(P_o)
    y = [0]*n
    for j in range(n/2):
        y[j] = y_e[j]+w^j*y_o[j]
        y[j+n/2] = y_e[j]-w^j*y_o[j]
    return y