FFT에서는 실수 자료형을 사용하기 때문에 실수 오차가 발생할 수 있고, 이는 즐거운 PS생활에 큰 지장을 줄 수 있습니다. 특히 FFT 문제에서 수가 너무 크기 때문에 M으로 나눈 나머지를 출력한다.와 유사한 문장이 나온다면 실수 오차를 만나는 것을 각오해야 합니다.
이 글에서는 M으로 나눈 나머지를 출력한다.같은 문장이 나오는 문제에서 약간의 속도를 희생해 FFT의 정확도를 높이는 방법과 정수만 이용해서 FFT를 하는 방법에 대해 다룹니다.
※ 이 글은 독자가 분할 정복을 이용한 FFT 알고리즘(쿨리-튜키 알고리즘)을 알고있다는 전제 하에 설명합니다.
쿨리-튜키 알고리즘을 이용해 다항식 곱셈을 하는 방법은 여기에서 볼 수 있습니다.
쿨리-튜키 알고리즘을 이용한 다항식 곱셈의 구현 예시는 여기에서 볼 수 있습니다.
정확도 높은 FFT
몇몇 문제에서는 두 개의 다항식을 곱한 결과가 너무 크기 때문에 어떤 수 $M$으로 나눈 나머지를 요구하는 문제들이 있습니다.
일반적인 FFT는 실수 자료형(complex<double>)을 이용하기 때문에 수가 커지면 실수 오차가 생길 수 밖에 없고, 모듈러를 취하는 대부분의 문제들은 실수 오차 이슈가 발생합니다. 이 단락에서는 다항식을 곱한 뒤 모듈러를 취할 때 흔히 발생하는 실수 오차를 피하는 방법을 다룹니다.
실수 오차가 생기는 이유는 실수 자료형이 정확하게 표현할 수 없을 정도로 수가 너무 커지기 때문입니다. 그러므로 계산 과정에서 나오는 수들의 크기가 작아지도록 조절해주면 실수 오차를 피할 수 있을 것입니다.
다항식 $A(x)$와 $B(x)$는 각각 아래와 같이 작은 다항식 2개로 쪼개줄 수 있습니다.
\[A(x) = A_1(x) + A_2(x) \times C \\ B(x) = B_1(x) + B_2(x) \times C\]두 n차 다항식 $A(x), B(x)$를 곱하는 것은 $A(x)B(x) = A_1(x)B_1(x)+C\times(A_1(x)B_2(x)+A_2(x)B_1(x))+C^2\times(A_2(x)B(x))$로 나타낼 수 있습니다.
$C \approx \sqrt M$ 정도로 정해주면 $A_1(x), A_2(x), B_1(x), B_2(x)$의 계수들의 최댓값이 $\sqrt M$ 정도이기 때문에 곱셈을 해도 계수들의 최댓값이 $Mn$ 정도가 되고, 기존의 방법보다 훨씬 작은 수들을 이용해 계산을 하기 때문에 실수 오차를 피할 수 있습니다.
여기 있는 코드는 $C = 2^{15}$일 때 위 알고리즘을 구현한 것입니다. 4번의 dft와 3번의 idft를 사용합니다.
koosaga님 github에 있는 다항식 라이브러리(링크)의 multiply_mod를 보면 2번의 dft와 2번의 idft만 이용해 하는 방법이 있는 것으로 보입니다.
정수 FFT
다항식 곱셈 결과를 특정 몇몇 소수로 나눈 나머지를 구할 때는 NTT라는 또 다른 방법을 사용할 수 있습니다.
$n = 2^k$차 다항식의 FFT를 구할 때 우리는 $w^n = 1$인 복소수 $w = \cos(2\pi/n) + i\times \sin(2\pi/n)$를 사용합니다. 이런 $w$는 sin, cos 등 실수로 구성된 복소수를 사용하기 때문에 실수 오차의 위험이 따라올 수 밖에 없습니다.
여기에서 생각해볼 수 있는 것은, $w^n ≡ 1\mod p$이면서 $w^0, w^1, \dots w^{n-1}$이 모두 서로 다른 $w$(원시근)와 $p$(소수)가 있고 이때 n이 $2^k$꼴의 배수라면, 원시근 $w$가 FFT에서 사용하는 복소수 $w$와 동일한 역할을 할 수 있게 됩니다.
결국, 적당히 큰 $2^k$꼴의 배수인 $n$과 조건을 만족하는 $w, p$를 찾을 수 있다면 $A(x)\times B(x)\mod p$를 정수만 이용해서 계산할 수 있습니다.
페르마의 소정리에 따르면 $w^{p-1} = 1\mod p$입니다. $n = p-1$은 $2^k$꼴의 배수가 되어야 하기 때문에, $p = a\times2^b +1$인 소수 $p$를 잡을 것입니다. 이런 $p$의 원시근 $w$를 찾았다면 $\cos(2\pi/n)+i\times\sin(2\pi/n)$ 대신 $w^{(p-1)/n}$을 이용해 FFT를 돌려줄 수 있습니다. $n ≤ 2^b$인 크기 $n$짜리 다항식에 대해 FFT mod P를 수행할 수 있습니다.
아래 표는 NTT에서 자주 사용하는 유명한 소수와 원시근입니다.
| p | a | b | w |
|---|---|---|---|
| 998,244,353 | 119 | 23 | 3 |
| 2,281,701,377 | 17 | 27 | 3 |
| 2,483,027,969 | 37 | 26 | 3 |
| 2,113,929,217 | 63 | 25 | 5 |
| 104,857,601 | 25 | 22 | 3 |
| 1,092,616,193 | 521 | 21 | 3 |
998,244,353과 비슷하게 생긴 993,244,853은 $248311213 \times 2^2 + 1$로, NTT를 하기에는 좋지 않은 수입니다. 문제에 가끔씩 나오기 때문에 주의해야 합니다.
여기 있는 코드는 10진수 자연수 곱셈을 NTT를 이용해 구현한 코드입니다.
NTT의 다양한 활용
(FFT의 결과) mod (square free number)
NTT는 $p = a\times 2^b + 1$이고 원시근 $w$가 존재할 때 (FFT의 결과) mod P를 정수만 이용해 계산할 수 있게 해줍니다. NTT와 중국인의 나머지 정리를 이용하면 (FFT의 결과) mod (square free number)도 해줄 수 있습니다.
$x \mod p_1 = a_1 \ x \mod p_2 = a_2 \ x \mod p_3 = a_3$
$p$가 square free number이고 이런 $(p_i, a_i)$쌍들이 주어졌을 때, $x$를 구할 수 있다면 $x \mod p$도 구할 수 있을 것입니다. 그리고 $x$는 중국인의 나머지 정리를 이용해 구할 수 있습니다.
NTT를 할 수 있는 여러 가지 소수로 NTT를 한 다음, 중국인의 나머지 정리를 이용해 결과를 합쳐주면 (FFT의 결과) mod (square free number)를 할 수 있습니다.
자연수 계수, 계수 < (1 « 63)
다항식의 계수가 자연수이고 등장하는 계수들이 long long 범위 안에 있으면서 mod를 하지 않는 경우, 다항식 곱셈을 NTT를 이용해서 할 수 있습니다.
계수가 $2^{63}$보다 작다는 것은 $\mod 2^{63}$을 해도 값이 달라지지 않는다는 것을 의미합니다. 그러므로 위에 나와있는 (FFT의 결과) mod (임의의 자연수)를 그대로 사용하면 됩니다.
10억~20억 근처의 소수 2개 정도 사용하면 답을 구할 수 있습니다.
C++로 구현한 코드는 여기에서 볼 수 있습니다.
참고 자료
- Codeforces Blog - General ideas (17. Modulo product of two polynomial with real-valued FFT)
- CP-Algorithms - FFT (FFT - Multiplication with arbitrary modulus)
- koosaga Blog - Fast Fourier Transform
- cubelover Blog - Number Theoretic Fast Fourier Transform