Prerequisite Knowledge

Polynomials

An expression of the form $$$P(x) = a_0 + a_1 x + \cdots + a_n x^n$$$ is called a polynomial.

The coefficient representation of a polynomial

We use the coefficient $$$a_0, a_1, \cdots, a_n$$$ to express the polynomial.

This is the most common form of polynomials in practice.

The point-value representation of a polynomial

Select $$$n + 1$$$ distinct values $$$x_0, x_1, \cdots, x_n$$$ and substitute them into the polynomial and calculate the corresponding values $$$P(x_0), P(x_1), \cdots, P(x_n)$$$, the set of all pairs

is called the point-value representation of the polynomial.

Complex Numbers

The imaginary unit $$$i$$$ is defined as the square root of $$$-1$$$, i.e. $$$i = \sqrt{-1}$$$.

A number of the form $$$a + bi (a, b \in \mathbb{R})$$$ is called a complex number.

Addition of Complex Numbers

$$$(a + bi) + (c + di) = (a + c) + (b + d)i$$$

Multiplication of Complex Numbers

$$$(a + bi)(c + di) = ac + adi + bci + bdi^2 = (ac - bd) + (ad + bc)i$$$

Another form of Complex Numbers

We can also express complex numbers in the form $$$r(\cos \theta + i \sin \theta)$$$ where $$$r$$$ is the distance from the origin to the point $$$(a, b)$$$ and $$$\theta$$$ is the angle formed by the vector $$$\overrightarrow{OZ} = (a, b)$$$ and the positive x-axis direction.

Roots of Unity

The $$$n$$$-th root of unity is defined as a complex number $$$z$$$ that satisfies the equation $$$z^n = 1$$$, where $$$n$$$ is a positive integer.

From the multiplication of complex numbers and trigonometric identities, we can find that $$$r_1(\cos \theta_1 + i \sin \theta_1) \cdot r_2(\cos \theta_2 + i \sin \theta_2) = r_1 r_2(\cos(\theta_1 + \theta_2) + i \sin (\theta_1 + \theta_2))$$$.

This then yields that $$$r^n (\cos(n \theta) + i \sin (n \theta)) = 1$$$ for all $$$n$$$-th root of unity.

So $$$r^n \cos(n \theta) = 1$$$ and $$$r^n \sin(n \theta) = 0$$$.

Squaring both sides gives $$$r^{2n} \cos^2(n \theta) = 1$$$ and $$$r^{2n} \sin^2(n \theta) = 0$$$.

Adding the two equations yields $$$r^{2n} (\cos^2(n \theta) + \sin^2(n \theta)) = r^{2n} = 1$$$.

Since, $$$r$$$ is a real number and $$$r \gt 0$$$, $$$r$$$ can be only $$$1$$$.

Which leads to $$$\cos(n \theta) = 1$$$ and $$$\sin(n \theta) = 0$$$.

So, $$$n \theta$$$ can only be $$$2 \pi k$$$ (Where $$$k$$$ is an integer).

Which means that $$$\displaystyle\theta = \frac{2 \pi k}{n}$$$.

So, all $$$n$$$-th roots of unity must lie on the unit circle and divide it into $$$n$$$ equal arcs.

In general, all $$$n$$$-th roots of unity are denoted by $$${\omega}_n^0, {\omega}_n^1, \cdots, {\omega}_n^{n-1}$$$.

From the argument above, we can conclude that $$$\displaystyle {\omega}_n^k = \cos \frac{2\pi k}{n} + i \sin \frac{2 \pi k}{n}$$$

Roots of unity have the following properties:

• $$$\displaystyle {\omega}_n^k = \cos \frac{2 \pi k}{n} + i \sin \frac{2 \pi k}{n} = \cos \frac{2 \pi \cdot 2k}{2n} + i \sin \frac{2 \pi \cdot 2k}{2n} = {\omega}_{2n}^{2k}$$$
• $$$\displaystyle {\omega}_{2n}^{n+k} = \cos \frac{2 \pi (n + k)}{2n} + i \sin \frac{2 \pi (n + k)}{2n} = \cos(\frac{2 \pi k}{2n} + \pi) + i \sin(\frac{2 \pi k}{2n} + \pi) = - \cos \frac{2 \pi k}{2n} - i \sin \frac{2 \pi k}{2n} = -{\omega}_{2n}^{k}$$$

Fast Fourier Transform (FFT)

The product (convolution) of two $$$n$$$-th degree polynomials is defined as:

For two $$$n$$$-th degree polynomials $$$F(x)$$$ and $$$G(x)$$$, the time complexity of naive multiplication using their coefficient representation is $$$\mathcal{O}(n^2)$$$.

But if we are using the point-value representation, the time complexity of is $$$\mathcal{O}(n)$$$ time.

Why?

Let $$$F(t)$$$ and $$$G(t)$$$ be the evaluations of $$$F(x)$$$ and $$$G(x)$$$ at $$$x = t$$$, then $$$(F * G)(t) = F(t)G(t)$$$ which can be computed in $$$\mathcal{O}(1)$$$ time.

Because there are $$$n + 1$$$ evaluations, so the total time complexity of multiplying two polynomials in point-value representation is $$$\mathcal{O}(n)$$$.

So if we find an efficient way of polynomial evaluation and interpolation we can optimize the complexity of doing convolutions.

Discrete Fourier Transform (DFT)

Consider an 8-th degree polynomial:

$$$F(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4 +a_5 x^5 + a_6 x^6 + a_7 x^7$$$

Let see what happens when we split $$$F(x)$$$ into even and odd degree terms:

$$$F(x) = (a_0 + a_2 x^2 + a_4 x^4 + a_6 x^6) + x(a_1 + a_3 x^2 + a_5 x^4 + a_7 x^6)$$$

Let's define two new polynomials:

$$$G(x) = a_0 + a_2 x + a_4 x^2 + a_6 x^3$$$

$$$H(x) = a_1 + a_3 x + a_5 x^2 + a_7x^3$$$

Then we have $$$F(x) = G(x^2) + x H(x^2)$$$.

Consider evaluating $$$F(x)$$$ at $$$x = \omega_n^k$$$ for all $$$k$$$ from $$$0$$$ to $$$n - 1$$$.

Since $$$\omega_{2n}^{2k} = \omega_n^k$$$, substituting $$$x = \omega_n^k$$$ for $$$\displaystyle k \lt \frac{n}{2}$$$ gives:

$$$\displaystyle F(\omega_n^k) = G(\omega_{n / 2}^k) + \omega_n^k H(\omega_{n / 2}^k)$$$

With the property $$$\omega_{2n}^{n + k} = -\omega_{2n}^k$$$ (the negative sign disappears after squaring), we get:

$$$\displaystyle F(\omega_n^{k + n / 2}) = G(\omega_{n / 2}^k) - \omega_n^k H(\omega_{n / 2}^k)$$$

Computing the point-value representation of $$$F(x)$$$ only requires computing the point-value representation of $$$G(x)$$$ and $$$H(x)$$$ (with problem size halved) and an $$$\mathcal{O}(n)$$$ merge step, leading to the recurrence:

$$$\displaystyle T(n) = 2T(\frac{n}{2}) + \mathcal{O}(n)$$$

By the master theorem, the time complexity of DFT is thus $$$\mathcal{O}(n\log n)$$$.

Inverse Discrete Fourier Transform (IDFT)

By using DFT, we can convert polynomials to point-value representations in $$$\mathcal{O}(n\log n)$$$ time and compute the point-value representation of their product in $$$\mathcal{O}(n)$$$ time.

Now the key problem is to interpolate the coefficient representation from the point-value representation.

We have to find $$$a_0, a_1, \cdots, a_{n - 1}$$$, known that $$$y_i = f(\omega_n^i)$$$.

Treat $$$y_i$$$ as the coefficients of a new polynomial:

$$$\displaystyle A(x) = \sum_{i = 0}^{n - 1} y_i x^i$$$

Consider computing the point-value representation of $$$A(x)$$$.

Let $$$b_i = w_n^{-i}$$$.

Define the sum $$$\displaystyle S(\omega_n^a) = \sum_{i = 0}^{n - 1}(\omega_n^a)^i$$$.

If $$$a = 0$$$, $$$S(\omega_n^a) = n$$$.

If $$$a \neq 0$$$, $$$S$$$ is the sum of a geometric sequence.

So, multiplying it by $$$\omega_n^a$$$, yields $$$\displaystyle \omega_n^a S(\omega_n^a) = \sum_{i = 1}^n (\omega_n^a)^i$$$.

Then subtracting the original sum, yields $$$(\omega_n^a - 1) S(\omega_n^a) = (\omega_n^a)^n - (\omega_n^a)^0 = 1 - 1 = 0$$$.

Substituting this back into the evaluation of $$$A(x)$$$ at $$$\omega_n^{-k}$$$, we can see that $$$A(\omega_n^{-k}) = \sum_{j=0}^{n-1}a_j S(\omega_n^{j-k})=n a_k$$$.

So, we take the root of unity as its reciprocal, and perform DFT on $$${y_0, y_1, \cdots, y_n}$$$, then dividing the output by $$$n$$$, we will obtain the coefficient representation of the product. Thus, we can do IDFT in $$$\mathcal{O}(n\log n)$$$ time.

Some Tips while doing FFT

For implementation convenience, resizing the polynomial to a degree that is a power of two (padding zeros for the coefficients beyond the original degree).

Precomputing the base-case index order

Let's simulate the spliting for an 8-th degree polynomial:

$$$x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7$$$

$$$x_0, x_2, x_4, x_6, x_1, x_3, x_5, x_7$$$

$$$x_0, x_4, x_2, x_6, x_1, x_5, x_3, x_7$$$

Writing indices in binary gives $$$000, 100, 010, 110, 001, 101, 011, 111$$$.

If you reverse the bits of the binary yields $$$000, 001, 010, 011, 100, 101, 110, 111$$$.

Precomputing this bit-reversal permutation eliminates the need to extract even/odd terms during divide-and-conquer.

No extra array allocation

For the merge step, store $$$\displaystyle G(\omega_{n / 2}^k)$$$ at the $$$k$$$-th position and $$$\displaystyle H(\omega_{n / 2}^k)$$$ at the $$$\displaystyle (k + \frac{n}{2})$$$-th position of the original array. This allows overwriting the array for $$$F(x)$$$'s point-values.

Number Theoretic Transform (NTT)

FFT has precision issues for large numbers and cannot handle modulo operations (which require integer arithmetic). NTT replaces complex roots of unity with primitive roots (because they have similar mathematical properties).

For a prime $$$p = qn + 1$$$ ($$$n = 2^m$$$), a primitive root $$$g$$$ of $$$p$$$ satisfies $$$g^{qn} \equiv 1 \pmod{p}$$$.

Define $$$g_n = g^q$$$, and all we need to do is to replace $$${\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}}$$$ with $$${g_n^0, g_n^1, \cdots, g_n^{n-1}}$$$ in all computations.

Common moduli for NTT with primitive root $$$g = 3$$$: $$$p = 998244353 = 7 \times 17 \times 2^{23} + 1, p = 1004535809 = 479 \times 2^{21} + 1$$$

Polynomial Differentiation

Given a polynomial $$$\displaystyle F(x) = \sum_{i = 0}^n a_i x^i$$$, its derivative is $$$\displaystyle F'(x) = \sum_{i = 0}^{n - 1} a_{i + 1}(i + 1) x^i$$$

Polynomial Integration

Given a polynomial $$$\displaystyle F(x) = \sum_{i = 0}^n a_i x^i$$$, its indefinite integral (with constant term $$$0$$$) is $$$\displaystyle\int F(x) dx=\sum_{i = 1}^{n + 1} \frac{a_{i - 1}}{i} x^i$$$.

Polynomial Inversion

The inverse of an $$$n$$$-th degree polynomial $$$F(x)$$$ is a polynomial $$$G(x)$$$ such that $$$F(x)G(x) \equiv 1 \pmod{x^n}$$$.

This can be solved via divide and conquer.

Suppose we have found the inverse $$$H(x)$$$ of $$$F(x)$$$ modulo $$$x^{n / 2}$$$, by definition $$$F(x)H(x) \equiv 1 \pmod{x^{n / 2}}$$$.

Since $$$F(x)G(x) \equiv 1 \pmod{x^{n / 2}}$$$, we have $$$G(x) - H(x) \equiv 0 \pmod{x^{n / 2}}$$$.

Squaring both sides gives $$$G^2(x) - 2G(x)H(x) + H^2(x) \equiv 0 \pmod{x^n}$$$.

Multiplying both sides by $$$F(x)$$$ yields $$$F(x)G^2(x) - 2F(x)G(x)H(x) + F(x)H^2(x) \equiv 0 \pmod{x^n}$$$.

Since $$$F(x)G(x) \equiv 1 \pmod{x^{n}}$$$, the equation simplifies to $$$G(x) - 2H(x) + F(x)H^2(x) \equiv 0 \pmod{x^n}$$$, so $$$G(x) \equiv 2H(x) - F(x)H^2(x) \pmod{x^n}$$$.

The base case of the process is a constant polynomial $$$c$$$: its inverse is the modular inverse of $$$c$$$.

Polynomial Logarithm

We are to find a polynomial $$$G(x)$$$ such that $$$G(x) \equiv \ln F(x) \pmod{x^n}$$$, where $$$\ln F(x)$$$ is defined by substituting $$$F(x)$$$ into the Taylor series of $$$\ln(1 + x)$$$.

By the chain rule of differentiation: $$$\displaystyle G'(x) = \frac{F'(x)}{F(x)}$$$.

So we can compute $$$G(x)$$$ as follows:

1. Compute the derivative $$$F'(x)$$$ of $$$F(x)$$$
2. Multiply $$$F'(x)$$$ by the inverse of $$$F(x)$$$ modulo $$$x^n$$$
3. Compute the indefinite integral of the product modulo $$$x^n$$$

Polynomial Exponential

We are to find a polynomial $$$G(x)$$$ such that $$$G(x) \equiv e^{F(x)} \pmod{x^n}$$$, where $$$e^{F(x)}$$$ is defined by substituting $$$F(x)$$$ into the Taylor series of $$$e^x$$$.

This can be solved via divide and conquer.

Using Newton's method, suppose we have found $$$H(x)$$$ as the exponential of $$$F(x)$$$ modulo $$$x^{n / 2}$$$ then the exponential modulo $$$x^n$$$ is given by $$$G(x) \equiv H(x) \bigl(1 - \ln H(x) + F(x) \bigr) \pmod{x^n}$$$,

Polynomial Powers

We are to find a polynomial $$$G(x)$$$ such that $$$G(x) \equiv F^m(x) \pmod{x^n}$$$ for a constant exponent $$$m$$$.

Taking the polynomial logarithm of both sides yields $$$\ln G(x) \equiv m \ln F(x) \pmod{x^n}$$$.

So we can compute $$$G(x)$$$ as follows:

1. Compute $$$\ln F(x)$$$ modulo $$$x^n$$$
2. Multiply the result by the constant $$$m$$$
3. Compute the polynomial exponential of the product modulo $$$x^n$$$

Some points that you should keep in mind:

• Polynomial logarithm requires the constant term of $$$F(x)$$$ to be $$$1$$$. If the constant term is $$$c \neq 1$$$, multiply all coefficients of $$$F(x)$$$ by the modular inverse of $$$c$$$
• If $$$F(x)$$$ has no constant term (i.e., $$$F(x) = x^c \cdot F_0(x)$$$ where $$$F_0(x)$$$ has a non-zero constant term), divide $$$F(x)$$$ by $$$x^c$$$, compute $$$F_0^m(x)$$$ modulo $$$x^{n - cm}$$$, then multiply the result by $$$x^{cm}$$$

Polynomial Square Root

We are to find a polynomial $$$G(x)$$$ such that $$$G^2(x) \equiv F(x) \pmod{x^n}$$$.

This can be solved via divide and conquer.

Suppose we have found the $$$H(x)$$$ as the square root of $$$F(x)$$$ modulo $$$x^{n / 2}$$$, by definition $$$H^2(x) \equiv F(x) \pmod{x^{n / 2}}$$$.

Since $$$G^2(x) \equiv F(x) \pmod{x^{n / 2}}$$$, we have $$$G(x) - H(x) \equiv 0 \pmod{x^{n / 2}}$$$.

Squaring both sides gives $$$G^2(x) - 2G(x)H(x) + H^2(x) \equiv 0 \pmod{x^n}$$$.

Substituting $$$G^2(x) \equiv F(x) \pmod{x^n}$$$ gives $$$F(x) - 2G(x)H(x) + H^2(x) \equiv 0 \pmod{x^n}$$$.

So, $$$\displaystyle G(x) \equiv \frac{H^2(x) + F(x)}{2H(x)} \pmod{x^n}$$$.

Polynomial Modular Reduction

Given an $$$n$$$-th degree polynomial $$$F(x)$$$ and an $$$m$$$-th degree polynomial $$$G(x) (n \geq m)$$$, find a quotient polynomial $$$Q(x)$$$ (degree $$$n - m$$$) and a remainder polynomial $$$R(x)$$$ (degree $$$ \lt m$$$) such that $$$F(x) = G(x)Q(x) + R(x)$$$.

Define the reverse polynomial of a polynomial $$$\displaystyle P(x) = \sum_{i = 0}^k a_i x^i$$$ as $$$\displaystyle P_r(x) = x^k P(\frac{1}{x}) = \sum_{i = 0}^k a_{k - i} x^i$$$ (its coefficients are the reverse of $$$P(x)$$$'s coefficients).

Substitute $$$x$$$ for $$$\displaystyle \frac{1}{x}$$$ into the division formula and multiply by $$$x^n$$$: $$$\displaystyle x^n F(\frac{1}{x}) = x^m G(\frac{1}{x}) \cdot x^{n - m} Q(\frac{1}{x}) + x^n R(\frac{1}{x})$$$.

Since $$$\displaystyle x^n R(\frac{1}{x}) = x^{n - m + 1} \cdot x^{m - 1} R(\frac{1}{x})$$$ is divisible by $$$x^{n - m + 1}$$$, the equation simplifies to: $$$F_r(x) \equiv G_r(x) \cdot Q_r(x) \pmod{x^{n - m + 1}}$$$.

The reverse quotient polynomial is thus $$$Q_r(x) \equiv F_r(x) \cdot G_r^{-1}(x) \pmod{x^{n - m + 1}}$$$.

So we can compute $$$Q(x)$$$ and $$$R(x)$$$ as follows:

1. Compute $$$F_r(x)$$$ and $$$G_r(x)$$$
2. Find the inverse of $$$G_r(x)$$$ modulo $$$x^{n - m + 1}$$$
3. Multiply to get $$$Q_r(x)$$$, then reverse its coefficients to get $$$Q(x)$$$
4. Compute the remainder $$$R(x) = F(x) - G(x)Q(x)$$$