<자료구조 알고리즘> FFT를 이용한 다항식 곱셈 알고리즘

 

다항식 곱셈은 수학에서 아주 기본적으로 하는 연산 중에 하나이다.
$$A(x) = 5x+2$$

$$B(x) = x^2+2x+3$$
이 두식을 곱한 C(x)를 구한다고 하면
$$C(x) = (5x+2)(x^2+2x+3) \\ = 5x \times x^2 + 5x \times 2x+2\times x^2 + 5x \times 3+2 \times 2x + 2 \times 3 = 5x^3 + 14x^2 + 19x + 6$$
이렇게 쉽게 계산할 수 있다. 하지만 컴퓨터로 짜려면 어떻게 해야 할까?

다항식을 나타내는 방법에는 대표적으로 coefficient representation과 value representation 이 있다.
coefficient representation으로 위의 A(x)와 B(x)를 나타내 보자면
$$A(x) = [5,2]$$
$$B(x) = [1,2,3]$$
으로 나타낼 수 있다.이를 간단하게 곱하면 $$C(x) = [5,14,19,6]$$을 구할 수 있다. 하지만 이러한 방식의 시간 복잡도는$O(n^2)$이 된다.

이를 조금 더 빠르게 할 수 있을까?

가장 중요한 것은 좌표평면 위에 n+1개의 점이 주어진다면 이 점들을 모두 지나는 n차 함수는 유일하게 정해진다는 것이다.
(자세한 증명은 여기서:https://math.stackexchange.com/questions/426932/why-are-vandermonde-matrices-invertible)

이를 이용하여 대략적인 방식을 간단하게 정리해 보자면

a-1차 다항식 A(x) 와 b-1차 다항식 B(x)를 곱하여 n-1차 다항식 C(x)를 구한다고 하자 (n = a + b - 1)
1. A(x)에 n개의 값을 대입하여 n개의 함숫값을 계산하고 B(x)에 n개의 값을 대입하여 n개의 함숫값을 계산하여 두 다항식을 value representaition으로 나타낸다.
2. 계산한 두 함수의 함수값을 곱하여 C(x)의 함수값을 구한다.
3. 구한 C(x)의 함숫값들을 이용하여 C(x)를 구한다.

 

 

1. FFT (A(x)에 n개의 값을 대입하여 n개의 함숫값을 계산하고 B(x)에 n개의 값을 대입하여 n개의 함숫값을 계산하여 두 다항식을 value representaition으로 나타낸다.)

첫번째 과정부터 해보도록 하겠다.
a차의 A(x)에 b차의 B(x)에 n개의 값을 단순히 대입하여 계산한다고 하면 $O(n^2)$의 시간복잡도가 걸린다. 이렇게 된다면 위의 방식에 비해 빠르지 않을 것이다. 그러면 어떻게 더 빠르게 할 수 있을까?

어떤 다항식에 한개의 값을 대입했을 때 바로 알 수 있는 다른 한 값이 무엇인지 생각해 보자
$x^8+x^6+x^4+x^2$ 이라는 식에 m을 대입하면 $m^8+m^6+m^4+m^2$된다 여기서 우리는 이 식에 $-m$을 대입했을 때의 값도 $m^8+m^6+m^4+m^2$일 것이라고 계산하지 않고도 알 수 있다.
x$^7+x^5+x^3+x^1$이라는 식에 n을 대입하면 $m^7+m^5+m^3+m^1$이 된다. 여기서도 우리는 이식에 $-m$ 을 대입했을 때의 값이 $-m^7-m^5-m^3-m^1$ 이라는 것을 -1만 곱하여 알 수 있다.
이러한 성질을 이용하면 다항식 A(x) B(x)에 n개의 값을 대입해 함숫값을 구하는 과정을 더 빠르게 할 수 있을 것이다.

간단하게 하기 위해 A(x)의 계수가 모두 1이라고 하겠다.
$$A(x) = x^a + x^{a-1} + x^{a-2} + \dots +x^2 + x + 1$$
위 식을 차수가 홀수인 항과 짝수인 항으로 나누어 보자
$$A(x) = (x^{a-1} + x^{a-3} + \dots + x^2+1 )+x(x^{a-1}+x^{a-3}+ \dots +x^2+1)$$
$$A(x) = P_0(x^2) + xP_1(x^2)$$
그러면 A(x)를 두 함수로 나타낼 수 있다.
이때 $$x_1, -x_1, x_2, -x_2, \dots  x_{a+1/2}, -x_{a+1/2}$$의 함숫값을 구한다고 한다면
$$A(x_1) = P_0(x_1^2) + x_1P_1(x_1^2)$$$$A(-x_1) = P_0(x_1^2) - x_1P_1(x_1^2)$$

$$A(x_2) = P_0(x_2^2) + x_2P_1(x_2^2)$$$$A(-x_2) = P_0(x_2^2) - x_2P_1(x_2^2)$$
$$.$$
$$.$$
$$.$$
$$A(x_{(a+1)/2}) = P_0(x_{(a+1)/2}^2) + x_{(a+1)/2}P_1(x_{(a+1)/2}^2)$$$$A(-x_{(a+1)/2}) = P_0(x_{(a+1)/2}^2) - x_{(a+1)/2}P_1(x_{(a+1)/2}^2)$$
으로 구할 수 있다. 여기서 $P_0(x_1^2)$의 값과 $P_1(x_1^2)$의 값을 다시 쓸 수 있는 것이다!!

그리고 이는 $P_0(x)$와 $P_1(x)$에 ${x_1}^2, {x_2}^2, \dots {x_{a+1}}^2$를 대입한 함숫값을 찾는 문제로 이어지고
다시 위의 방법을 사용하여 재귀적으로 구할 수 있을 것 같다.
그런데 여기서 문제가 생겼다.
위에서 속도가 빨라질 수 있었던 이유는 +-짝을 이룬 값들을 이용하였기 때문인데 ${x_1}^2, {x_2}^2, \dots {x_{a+1}}^2$도 +-짝으로 이루어 질수 있을까?


N-th root of uinity
를 이용하면 된다.

 

 

이 수들은 다음과 같은 성질을 가진다

즉 항상 원점 대칭인 점이 있다는 이야기 이고 이 모든 수를 제곱하여도 절댓값이 같고 부호가 다른 수들끼리 짝을 지을 수 있다는 이야기이다. 이 수들을 대입하여 함숫값을 구하면 위에 방법을 통해 재귀적으로 구할 수 있다. 

(단, n는 2의 제곱수)
$$x_k = \omega_n^k = e^\frac{2{\pi}ik}{n} = cos(\frac{2{\pi}k}{n}) + sin(\frac{2{\pi}k}{n})$$
$$\omega_n^k = -\omega_n^{k+\frac{n}{2}}$$

그런데 만약 n이 2의 제곱수가 아니라면 어떻게 해야 할까? 더 많은 값을 집어 넣는 것은 상관이 없기 때문에 n을 보다 큰 2의 제곱수로 만들고 함숫값을 구해주면 된다. 예를들어 3차 다항식 두개를 곱한다고 하면 곱했을 때 6차 다항식이 나오므로 n을 8이라고 하고 8개의 값으로 진행하면 된다. 이렇게 첫 번째 단계가 끝났다. 

 

시간복잡도 O(n log n)

2. 계산한 두 함수의 함수값을 곱하여 C(x)의 함수값을 구한다.

fft를 이용해 value representation 으로 바꾼 다항식 두개를 곱해준다. 

시간복잡도 O(n)

(단, n는 2의 제곱수)
$$x_k = \omega_n^k = e^\frac{2{\pi}ik}{n} = cos(\frac{2{\pi}k}{n}) + sin(\frac{2{\pi}k}{n})$$
$$\omega_n^k = -\omega_n^{k+\frac{n}{2}}$$

 

3. IFFT (구한 C(x)의 함숫값들을 이용하여 C(x)를 구한다.)

마지막 역변환 단계이다.

$$C(x) = \frac{1}{n}(f(\bar{x_{n-1}})x^{n-1} + f(\bar{x_{n-2}})x^{n-2} + \dots + f(\bar{x_1})x + f(\bar{x_0}))$$

으로 구할 수 있다.

시간복잡도 O(n)

 

어떤 함수 f(x)를 푸리에 변환하는 과정을 행렬로 표현해 보자면 다음과 같다.($C_k$는 k차 항의 계수)

 

 

$$\begin{pmatrix} 1 & 1 & 1 & 1 & \dots & 1 & 1 \\ 1 & \omega & \omega^2 & \omega^3 & \dots & \omega^{n-2} & \omega^{n-1} \\ 1 & \omega^2 & \omega^4 & \omega^6 & \dots & \omega^{2(n-2)} & \omega^{2(n-1)} \\ \vdots &  \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 1 & \omega^{(n-2)} & \omega^{2(n-2)} & \omega^{3(n-2)} & \dots & \omega^{(n-2)(n-2)} & \omega^{(n-1)(n-2)} & \end{pmatrix} \begin{bmatrix} C_0 \\ C_1 \\ C_2 \\ \vdots \\ C_{n-1} \end{bmatrix} = \begin{bmatrix} f(x_0) \\ f(x_1) \\ f(x_2) \\ \vdots \\ f(x_{n-1}) \end{bmatrix}$$

 

 

$$\begin{bmatrix} C_0 \\ C_1 \\ C_2 \\ \vdots \\ C_{n-1} \end{bmatrix} = \begin{pmatrix} 1 & 1 & 1 & 1 & \dots & 1 & 1 \\ 1 & \omega & \omega^2 & \omega^3 & \dots & \omega^{n-2} & \omega^{n-1} \\ 1 & \omega^2 & \omega^4 & \omega^6 & \dots & \omega^{2(n-2)} & \omega^{2(n-1)} \\ \vdots &  \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 1 & \omega^{(n-2)} & \omega^{2(n-2)} & \omega^{3(n-2)} & \dots & \omega^{(n-2)(n-2)} & \omega^{(n-1)(n-2)} & \end{pmatrix}^{-1} \begin{bmatrix} f(x_0) \\ f(x_1) \\ f(x_2) \\ \vdots \\ f(x_{n-1}) \end{bmatrix}$$

즉$$ V · \bar{V} = n · I$$

라는 것만 증명하면 되는 것이다.

$$P = V ·\bar{V}$$

이라고 하자

 $$P_ij = (\text{row i of }V)(\text{column j of }\bar{V})$$

 

(저 오메가들의 행렬을 V라고 하자)

$P_{ij} = (\text{row i of }V)(\text{column j of} \bar{V})$

$= \displaystyle\sum_{m=0}^{n-1} e^{ij2\pi \frac{m}{n}} \cdot e^{-i2\pi k\frac{m}{n}} = \displaystyle\sum_{m=0}^{n-1} e^{i2\pi \frac{m}{n}(j-k)}$

 

j = k 일 때,

$P_{ij} = n$

j != k 일때

$P_{ij} = \displaystyle\sum_{m=0}^{n-1} (e^{i2\pi \frac{(j-k)}{n}})^m$

$ = \frac{(e^{i2\pi \frac{(j-k)}{n}})^n - 1}{e^{i2\pi \frac{(j-k)}{n}} -1} = \frac{1-1}{e^{i2\pi \frac{(j-k)}{n}} -1} = 0$

따라서 $\bar{V}/n$은 V의 역행렬이다.

 

이로써 증명이 완료 되었다.

이 또한 같은 방식을 사용하므로 시간복잡도는 O(n log n)이다.