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

 

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

$$B(x)=x^2+2x+3$$
이 두식을 곱한 C(x)를 구한다고 하면
$$\begin{gathered} 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 \end{gathered}$$
이렇게 쉽게 계산할 수 있다. 하지만 컴퓨터로 짜려면 어떻게 해야 할까?

다항식을 나타내는 방법에는 대표적으로 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}+\cdots +x^2+x+1$$
위 식을 차수가 홀수인 항과 짝수인 항으로 나누어 보자
$$A(x)=(x^{a-1}+x^{a-3}+\cdots +x^2+1)+x(x^{a-1}+x^{a-3}+\cdots +x^2+1)$$
$$A(x)=P_0(x^2)+x{P}_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_1{P}_1(x_1^2)$$$$A(-x_1)=P_0(x_1^2)-x_1{P}_1(x_1^2)$$

$$A(x_2)=P_0(x_2^2)+x_2{P}_1(x_2^2)$$$$A(-x_2)=P_0(x_2^2)-x_2{P}_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(\overline{x_{n-1}})x^{n-1}+f(\overline{x_{n-2}})x^{n-2}+\cdots +f(\overline{x_1})x+f(\overline{x_0}))$$

으로 구할 수 있다.

시간복잡도 O(n)

 

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

 

 

$$\left(\begin{array}{ccccccc}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)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& {\omega }^{(n-2)}& {\omega }^{2(n-2)}& {\omega }^{3(n-2)}& \dots & {\omega }^{(n-2)(n-2)}& {\omega }^{(n-1)(n-2)}& \end{array}\right)\left[\begin{array}{c}{C}_0\\ C_1\\ C_2\\ ⋮\\ C_{n-1}\end{array}\right]=\left[\begin{array}{c}f(x_0)\\ f(x_1)\\ f(x_2)\\ ⋮\\ f(x_{n-1})\end{array}\right]$$

 

 

$$\left[\begin{array}{c}{C}_0\\ C_1\\ C_2\\ ⋮\\ C_{n-1}\end{array}\right]={\left(\begin{array}{ccccccc}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)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& {\omega }^{(n-2)}& {\omega }^{2(n-2)}& {\omega }^{3(n-2)}& \dots & {\omega }^{(n-2)(n-2)}& {\omega }^{(n-1)(n-2)}& \end{array}\right)}^{-1}\left[\begin{array}{c}f(x_0)\\ f(x_1)\\ f(x_2)\\ ⋮\\ f(x_{n-1})\end{array}\right]$$

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

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

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

이라고 하자

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

 

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

$P_{ij}=(\text{row i of }V)(\text{column j of}\overline{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$

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

 

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

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