FFT

최대 차수가 \(n\)인 두 다항식 \(A, B\)가 있습니다.

두 다항식의 곱 \(A \times B\)을 계산하는 작업은 일반적으로 \(O(n^2)\)의 시간복잡도를 가지게 됩니다.

이 작업을 무려 \(O(n\log n)\)에 할 수 있는 FFT 알고리즘에 대해 알아봅시다.

 

 

 

먼저, 다항식을 나타내는 2가지 방법에 대해 알아봅시다.

\(n\)차 다항식을 나타낸다고 생각해 봅시다.

먼저, 일반적으로 \(n\)차 다항식의 \(n+1\)개의 계수를 차례대로 늘어놓는 방법이 있습니다.

 

$$ax^2+bx+c \rightarrow \{c,b,a\}$$

\(i\)번째 원소는 \(i\)차 항을 나타내기 때문에 역순입니다.

이를 Coefficient Representation이라고 합니다.

 

두 번째로, 이 다항식에 \(n+1\)개의 서로 다른 \(x\)값을 대입해 \(f(x)\)를 계산하여, \(n+1\)개의 쌍으로 나타내는 방법이 있습니다.

$$ax^2+bx+c \rightarrow \{(x_0, f(x_0)),(x_1, f(x_1)),(x_2, f(x_2))\}$$

이를 Value Representation이라고 합니다.

 

이 두 번째 방법을 적극적으로 이용해 봅시다.

 

다항식 \(A\)와 \(B\)가 각각 \(n\)차 다항식이라고 하고, 이 둘의 곱을 \(C\)라고 가정합시다.

\(C\)는 \(2n\)차 다항식이므로, 만약 우리가 \(2n+1\)개의 서로 다른 \((x, C(x))\)쌍이 있으면, 이 값을 이용해 \(C\)를 복원할 수 있게 됩니다.

 

이 방법을 이용하여,

\(A\)와 \(B\)를 \(2n+1\)개의 서로 다른 \(x\)에 대해, \((x_i, A(x_i))\)과 \((x_i, B(x_i))\)쌍으로 바꾼 다음,

함수값을 곱하여 \(2n+1\)개의 \((x_i, C(x_i))\)쌍으로 만들고, 이를 \(C\)로 복원해 낼 수 있게 됩니다.

 

 

이를 위해 2가지를 알아야 합니다.

첫 번째는 어떤 Coefficient Representation으로 나타내어진 다항식을 Value Representation으로 변환하는 방법,

두 번째는 반대로 Value Representation으로 나타내어진 다항식을 Coefficient Representation으로 변환하는 방법입니다.

 

이 변환에 FFT가 쓰이게 됩니다.

 

 

먼저, Coefficient Representation을 Value Representation으로 변환하는 방법에 대해 알아봅시다.

우리는 \(2n+1\)개의 서로 다른 \(x\)에 대해, 함수값을 알아내야 합니다.

당연히 이를 나이브하게 구하려면 \(O(n^2)\)이므로, 조금 더 빠른 방법을 생각해 봅시다.

 

만약 우리가 \(A\)라는 다항식이 차수가 짝수인 항으로만 이루어져 있다고 가정해 봅시다.

ex) \(A(x) = x^4 + 3x^2 + 2\)

여기에서 우리가 어떤 \((x, A(x))\)쌍 하나를 이미 알고 있다면, 이를 이용해 \((-x, A(x))\)라는 쌍을 하나 더 알 수 있습니다.

 

반대로, \(B\)라는 다항식이 차수가 홀수인 항으로만 이루어져 있다고 가정해 봅시다.

ex) \(B(x) = 2x^5 + 3x^3 + 7x\)

역시 우리가 어떤 \((x, B(x))\)라는 쌍 하나를 알고 있다면, 이를 이용해 \((-x, -B(x))\)라는 쌍을 하나 더 알 수 있습니다.

 

그러면 일반적인 다항식을 짝수 차수항과 홀수 차수항 둘로 나누어서 각각 \(n\)개의 서로 다른 \(x\)에 대해 함수값을 알아낸다면,

이를 이용해 \(2n\)개의 함수값을 알 수 있다는 뜻이 됩니다.

 

정리하면, 어떤 일반적인 \(n\)차 다항식을 \(P\)라고 했을 때,

(\(P_e : P\)의 짝수 차수항만 모아놓은 다항식, \(P_o : P\)의 홀수 차수항만 모아놓은 다항식)

$$P(x) = P_e(x^2) + xP_o(x^2)$$

이므로, 우리가 어떤 \(x_i\)에 대해 \(P_e(x_i^2)\)와 \(P_o(x_i^2)\)를 알고 있으면,

$$P(x_i) = P_e(x_i^2) + x_iP_o(x_i^2)$$

$$P(-x_i) = P_e(x_i^2) - x_iP_o(x_i^2)$$

이와 같이 \(x_i\)와 \(-x_i\)에 대한 함수값도 알 수 있게 됩니다.

 

그런데 이를 당연하게도, \(P_e\)와 \(P_o\)에도 적용할 수 있습니다.

이를 재귀함수로 구현한다고 생각해 보면,

\(P_e\)와 \(P_o\)는 \(P\)에 비해 최고차항의 차수가 절반으로 줄어든 상태이고,

둘을 합치는 과정에서 덧셈의 횟수는 최고차항의 차수와 비례하므로,

우리는 이 알고리즘의 시간복잡도가 \(O(n\log n)\)이 될 것이라는 사실을 알 수 있습니다.

 

그러면 맨 처음 \(P\)를 Value Representation으로 바꾸기 위해 \(n\)개의 서로 다른 \(x\)값을 넣는다고 했는데, 이 \(x\)값이 어떤 값이어야 하는지 알아봅시다.

이 재귀 함수의 기저사례는 다항식의 차수가 0차일 때일 것입니다.

ex) \(Y(x) = p_0\)

(별로 의미없지만) \(x=1\)을 대입해서 \(Y(1) = p_0\)을 만든 다음, 이 값을 반환합니다.

이 다음 단계를 생각해 봅시다.

우리는 지금 \(x_i^2\)에 대한 함수값을 알아냈고, 이를 이용해 \(x_i\)과 \(-x_i\)의 함수값을 계산해야 합니다.

우리는 \(x_i^2 = 1\)에 대한 함수값을 알아냈으니, \(x_i = 1, x_i = -1\)에 대한 함수값을 알아낼 수 있습니다. 이 값들을 반환합니다.

그 다음 단계를 생각해봅시다.

1은 이전과 다를바가 없을 것인데, 문제는 -1입니다.

\(x_i^2 = -1\)에 대한 함수값을 알아냈을 때, \(x_i\)과 \(-x_i\)의 함수값을 알아내야 합니다.

-1의 두 제곱근은 무엇인가요? \(i\)와 \(-i\)입니다.

그렇습니다. 우리는 실수 다항식에 대한 계산을 하기 위해 \(x\)에 복소수를 대입해야 합니다. (!)

 

일반적으로, \(n\)이 \(2^k\) (\(k\)는 자연수)꼴일 때,

\(n-1\)차 다항식 \(P\)를 구하기 위해 대입해야 하는 \(x\)들의 값은 \(x^n = 1\)이라는 식의 해입니다. (복소수 해 포함)

 

이 해들을 구하는 방법은 이 글에서는 생략하고, 아무튼 그 값들은 다음과 같습니다.

\(\omega = \cos(\frac{2\pi}{n}) + i\sin(\frac{2\pi}{n})\)일 때,

$$\{\omega^0, \omega^1, \omega^2, \cdots, \omega^{n-1}\}$$

 

그러면 이 각각의 \(x\)값들의 함수값 \(P(x)\)의 값을 계산할 수 있습니다.

 

 

지금까지 우리는 Coefficient Representation을 Value Representation으로 변환하는 방법인 FFT에 대해 알아봤습니다.

그러면 반대로, Value Representation을 Coefficient Representation으로는 어떻게 변환할 수 있을까요?

 

우리가 지금까지 만들었던 Value Representation을 행렬로 나타내면 다음과 같습니다.

\(p_i\)는 \(P\)의 \(i\)차항의 계수입니다.

$$
\begin{pmatrix}
P(\omega^0) \\ P(\omega^1) \\ P(\omega^2) \\ \vdots \\ P(\omega^{n-1})
\end{pmatrix}=\begin{pmatrix}
\omega^0 & \omega^0 & \omega^0 & \dots & \omega^0 \\ 
\omega^0 & \omega^1 & \omega^2 & \dots & \omega^{n-1} \\ 
\omega^0 & \omega^2 & \omega^4 & \dots & \omega^{2(n-1)} \\ 
\vdots & \vdots & \vdots & \ddots & \vdots \\ 
\omega^0 & \omega^{n-1} & \omega^{2(n-1)} & \dots & \omega^{(n-1)(n-1)}
\end{pmatrix}
\begin{pmatrix}
p_0 \\ p_1 \\ p_2 \\ \vdots \\ p_{n-1}
\end{pmatrix}
$$

 

우리는 지금까지 \(p_i\)를 이용해 \(P(\omega^i)\)를 구했습니다.

반대로 \(P(\omega^i)\)를 이용해 \(p_i\)를 구하려면, 가운데에 있는 행렬의 역행렬을 구해 양변에 곱해주면 될 것입니다.

 

$$
\begin{pmatrix}
p_0 \\ p_1 \\ p_2 \\ \vdots \\ p_{n-1}
\end{pmatrix}=\begin{pmatrix}
\omega^0 & \omega^0 & \omega^0 & \dots & \omega^0 \\ 
\omega^0 & \omega^1 & \omega^2 & \dots & \omega^{n-1} \\ 
\omega^0 & \omega^2 & \omega^4 & \dots & \omega^{2(n-1)} \\ 
\vdots & \vdots & \vdots & \ddots & \vdots \\ 
\omega^0 & \omega^{n-1} & \omega^{2(n-1)} & \dots & \omega^{(n-1)(n-1)}
\end{pmatrix}^{-1}
\begin{pmatrix}
P(\omega^0) \\ P(\omega^1) \\ P(\omega^2) \\ \vdots \\ P(\omega^{n-1})
\end{pmatrix}
$$

 

이 가운데에 있는 행렬을 DFT행렬이라고 하는데, 이 행렬의 역행렬은 다음과 같음이 잘 알려져 있습니다.

 

$$
\begin{pmatrix}
\omega^0 & \omega^0 & \omega^0 & \dots & \omega^0 \\ 
\omega^0 & \omega^1 & \omega^2 & \dots & \omega^{n-1} \\ 
\omega^0 & \omega^2 & \omega^4 & \dots & \omega^{2(n-1)} \\ 
\vdots & \vdots & \vdots & \ddots & \vdots \\ 
\omega^0 & \omega^{n-1} & \omega^{2(n-1)} & \dots & \omega^{(n-1)(n-1)}
\end{pmatrix}^{-1}=\frac 1n
\begin{pmatrix}
\omega^0 & \omega^0 & \omega^0 & \dots & \omega^0 \\ 
\omega^0 & \omega^{-1} & \omega^{-2} & \dots & \omega^{-(n-1)} \\ 
\omega^0 & \omega^{-2} & \omega^{-4} & \dots & \omega^{-2(n-1)} \\ 
\vdots & \vdots & \vdots & \ddots & \vdots \\ 
\omega^0 & \omega^{-(n-1)} & \omega^{-2(n-1)} & \dots & \omega^{-(n-1)(n-1)}
\end{pmatrix}
$$

 

행렬의 각 원소의 \(\omega\)를 \(\omega^{-1}\)로 치환하고, 결과에 \(\frac 1n\)을 곱하면 사실상 같은 연산이라는 사실을 알 수 있습니다.

다시 말하면, 역 FFT는 FFT와 동일한 코드에서 \(\omega\)의 값만 조금 변경함으로써 구할 수 있다는 뜻입니다.

 

 

여기까지의 내용을 모두 이해하셨다면, 적어도 FFT를 코드로 짤 수는 있게 되었습니다!

 

www.acmicpc.net/problem/15576

15576번: 큰 수 곱셈 (2)

어떤 자연수 \(\overline{abcde}\)를 \(ax^4+bx^3+cx^2+dx+e\)라는 다항식에 \(x=10\)을 대입한 결과라고 생각할 수 있습니다.

따라서 FFT를 이용해 두 다항식을 곱한 다음 \(x\)에 값을 대입하는 방식으로,

\(n\)자리 수의 곱셈을 \(O(n\log n)\)에 처리할 수 있습니다.

 

#include <bits/stdc++.h>
using namespace std;
 
typedef long long ll;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef complex<double> cpdb;
 
ll gcd(ll a, ll b) { for (; b; a %= b, swap(a, b)); return a; }
 
const double PI = acos(-1);
string A, B;
 
cpdb w;
vector <cpdb> a, b;
 
void FFT(vector <cpdb>& p, cpdb w)
{
    int n = p.size();
    if (n == 1) return; // 기저 사례
 
    // P를 짝수 차수와 홀수 차수 항으로 분리
    vector <cpdb> even(n >> 1), odd(n >> 1);
    for (int i = 0; i < n; i++)
    {
        if (i & 1) odd[i >> 1] = p[i];
        else even[i >> 1] = p[i];
    }
 
    // 각각 FFT변환
    FFT(even, w * w), FFT(odd, w * w);
 
    cpdb cw(1, 0);
    for (int i = 0; i < n / 2; i++)
    {
        p[i] = even[i] + cw * odd[i];
        p[i + n / 2] = even[i] - cw * odd[i];
 
        cw *= w;
    }
}
 
vector <cpdb> mul(vector <cpdb>& a, vector <cpdb>& b)
{
    int n = 1;
    while (n < a.size() || n < b.size()) n <<= 1;
    n <<= 1;
 
    // FFT를 편하게 하기 위해 a, b의 크기를 2^k꼴로 바꿔준다
    a.resize(n), b.resize(n);
 
    cpdb w(cos(2 * PI / n), sin(2 * PI / n)); // w의 값 계산
 
    // a, b를 각각 FFT변환
    FFT(a, w); FFT(b, w);
 
    vector <cpdb> c(n);
    for (int i = 0; i < n; i++) c[i] = a[i] * b[i];
 
    // c를 FFT 역변환
    cpdb rw(cos(-2 * PI / n), sin(-2 * PI / n));
    FFT(c, rw);
    for (int i = 0; i < n; i++) c[i] /= n;
 
    return c;
}
 
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0), cout.tie(0);
 
    cin >> A >> B;
    for (int i = A.size() - 1; i >= 0; i--) a.push_back(A[i] - '0');
    for (int i = B.size() - 1; i >= 0; i--) b.push_back(B[i] - '0');
 
    vector <cpdb> c = mul(a, b);
 
    vector <int> res(c.size());
    for (int i = 0; i < c.size(); i++)
    {
        res[i] = (int)(c[i].real() + 0.5); // 실수오차가 있을 수 있으므로 반올림한다.
    }
 
    for (int i = 0; i + 1 < res.size(); i++)
    {
        res[i + 1] += res[i] / 10;
        res[i] %= 10;
    }
 
    while (!res.empty() && res.back() == 0) res.pop_back();
    if (res.empty()) res.push_back(0);
    reverse(res.begin(), res.end());
 
    for (int x : res) cout << x;
}

---

제 그룹의 문제집에서 연습 문제들을 관리하고 있습니다.
문제집의 문제들을 보고 싶으시다면, 가입 신청을 해 주세요.

 

www.acmicpc.net/group/7712

ANZ1217