Recently I've been thinking about alternatives to the usual Cooley-Tukey FFT for multiplying polynomials. I think I've come up with a pretty nifty algorithm for multiplying polynomials. So my question is, does anyone know if this FFT algorithm has name? Or is it new algorithm?
My algorithm calculates the product of two polynomials $P$ and $Q$ $\mod x^n - \omega$ in $O(n \log n)$ time, assuming $n$ is a power of two, and $\omega$ is a non-zero complex number.
Part 1. Main idea:
Suppose you know some polynomial $P \mod \, (x^{2n} - \omega)$, where $n$ is a non-negative integer, and $\omega$ is a complex number. Let $P_1 = P \mod (x^n + \sqrt{\omega})$ and $P_2 = P \mod (x^n - \sqrt{\omega})$. I think of $P_1$ and $P_2$ as two different "folds" of $P \mod (x^{2n} - \omega)$. Note that this folding operation is invertible. From knowing $P \mod (x^n + \sqrt{\omega})$ and $P \mod (x^n - \sqrt{\omega})$ you can "unfold" them into $P \mod (x^{2n} - \omega)$ by using the sum of the folds, and the difference of the folds.
import cmath
def fold(P, w):
nhalf = len(P)//2
return [P[i] + w * P[i + nhalf] for i in range(nhalf)]
def unfold(P1, P2, w):
return [(p1 + p2)/2 for p1, p2 in zip(P1, P2)] + \
[(p1 - p2)/(2 * w) for p1, p2 in zip(P1, P2)]
# Dummy identity function
def fold_and_unfold(P, w):
sqrtw = cmath.sqrt(w)
P1 = fold(P, sqrtw)
P2 = fold(P, -sqrtw)
return unfold(P1, P2, sqrtw)
print(fold_and_unfold([1,2,3,4], 1j + 100)) # This prints [1,2,3,4]
Part 2. Polynomial multiplication algorithm:
The key insight here is that if you want to efficiently multiply two polynomials $P$ and $Q$, then you can start by repeatedly folding $P$ and $Q$, then multiply their folded representations, and then unfold the result.
def multiply(P, Q, w):
"""
O(n log n) polynomial multiplication algorithm
Input:
Polynomial P and Q of length n, n must be a power of 2
w any non-zero complex number
Output:
Polynomial R = P * Q mod x^n - w
R has length n
"""
if len(P) == 1:
return [P[0] * Q[0]]
sqrtw = cmath.sqrt(w)
P1 = fold(P, sqrtw)
P2 = fold(P, -sqrtw)
Q1 = fold(Q, sqrtw)
Q2 = fold(Q, -sqrtw)
R1 = multiply(P1, Q1, sqrtw)
R2 = multiply(P2, Q2, -sqrtw)
return unfold(R1, R2, sqrtw)
# This prints[3, 10, 22, 40, 43, 38, 24]
print(multiply([1,2,3,4,0,0,0,0], [3,4,5,6,0,0,0,0], 1j + 100))
This algorithm is similar to Cooley-Tukey FFT based polynomial multiplication algorithms, and it uses exactly the same number of complex multiplications.
But there are some key difference. My algorithm works $\mod x^n - \omega$ for any non-zero complex $\omega$. Also, during the folding operations $\omega$ stays constant. This is very different to Cooley-Tukey FFT where nothing stays constant.