Motivation
Let’s take a look at how we multiply two polynomials of degree
We can use the distributive property to multiply two polynomials and then sum up coefficients for identical terms.
However, this leads to
Polynomial Representations
A polynomial of degree
- Coefficient Representation:
coefficients - Value Representation:
sample points
Note that any
Polynomial Multiplication
We observe that computing polynomial multiplication using the coefficients leads to
These are the polynomials f(x) and g(x) that we want to multiply:
Step1: Sample
Step2: Compute the pairwise product of samples.
Step3: Find the unique polynomial of degree
This new red curve is our polynomial
Now we need two magical operations to take us from the coefficient representation to the value representation and vice versa in less than
Fast Fourier Transform
Evaluation (DFT)
We know the result of multiplying two degree
Observation
If we want to sample the function
Arbitrary Polynomials
Let’s say we want to evaluate the function
Property 1: If we compute F(x) in this fashion, computing F(-x) becomes trivial.
Property 2: Note that
Recap: We start with a function
Challenge: The problem is as we go to recurse to the next step, when we decide to split our samples into positive and negative pairs
Complex Samples
We observed that if we pick a sample pair
Nth Roots of Unity: In general, the solutions to the equation
Basically, we are going to evaluate our polynomial at values
Implementation
We can implement this idea recursively in 10 lines of code. This implementation allows for evaluating
import math, numpy as np
def FFT(c): # coefficients: [-1, 2, 3, 0] for 3x^2+2x-1
n = len(c) # n must be a power of 2
if n == 1:
return c
w = np.exp(2j * math.pi / n)
fe, fo = FFT(c[::2]), FFT(c[1::2])
f = [0] * n
for i in range(n//2):
f[i] = fe[i] + w**i * fo[i]
f[i + n//2] = fe[i] - w**i * fo[i]
return f
Interpolation
Having transformed the coefficients into points and performed pairwise multiplication, the next step involves determining the polynomial that passes through these points, a process known as “interpolation”. Although interpolation may initially appear to be a more challenging task, we will explore the close relationship between evaluation and interpolation, leading to a surprisingly simple solution.
DFT Matrix
Note that while sampling converts a polynomial (coefficients) to a set of points, interpolation acts as the inverse of that function converting points to the polynomial. Since sampling is a linear transformation, naturally the next step would be write the sampling/evaluation operation as a matrix and try to find its inverse to achieve the interpolation matrix.
Let’s take another look at sampling/evaluation for polynomial
And in the FFT algorithm, the i-th evaluation point is the corresponding root of unity (
Inverse DFT Matrix
To be more precise, the DFT matrix also has a normalization factor of
This normalization factor is often added to achieve a unitary matrix. Considering that the inverse of a unitary matrix is the conjugate transpose, and given our matrix is symmetric, we can achieve this by computing the conjugate of each term (negating the imaginary part). For the specific case of
Yet, it’s crucial to recognize that our initial DFT transformation was essentially an evaluation, devoid of the
Multiplying this matrix by the evaluated points will yield the coefficients of the desired polynomial.
Now, comparing this with the original
- negative powers
- multiply by
at the end
Full Implementation
With a one line change in the original implementation of our FFT
function, we can have it compute inverse FFT as well. Just keep in mind that we have to divide the final result by IFFT
).
import math, numpy as np
def FFT(c, inv=False):
n = len(c) # coefficients: [-1, 2, 3, 0] for 3x^2+2x-1
if n == 1: # n must be a power of 2
return c
w = np.exp(2j * math.pi / n * (-1 if inv else 1))
fe, fo = FFT(c[::2], inv), FFT(c[1::2], inv)
f = [0] * n
for i in range(n//2):
f[i] = fe[i] + w**i * fo[i]
f[i + n//2] = fe[i] - w**i * fo[i]
return f
def IFFT(c):
f = FFT(c, inv=True)
return [i / len(c) for i in f]
Testing
Let’s put our original polynomials to the test:
Input
Note that the coefficients need to be sorted from the lowest power to the highest.
poly1 = [3, 4, -6, -2, 4]
poly2 = [6, -1, -9, 11, -1]
Code
Just keep in mind the output of IFFT is complex numbers, so we need to extract the real parts and round them to the nearest integer, since the types are float.
def make_power_of_two(arr):
target_len = 2 ** (len(arr) - 1).bit_length()
return arr + [0] * (target_len - len(arr))
def multiply_polynomials(poly1, poly2):
# make them equal length and powers of 2
target_len = 2 ** (len(poly1) + len(poly2) - 1).bit_length()
poly1 += [0] * (target_len - len(poly1))
poly2 += [0] * (target_len - len(poly2))
p1, p2 = FFT(poly1), FFT(poly2)
p_result = [p1[i] * p2[i] for i in range(target_len)]
return IFFT(p_result)
Output
[18, 21, -67, -9, 121, -56, -52, 46, -4, 0, 0, 0, 0, 0, 0, 0]