모듈러 «나눗셈» - 알고리즘카

올림피아드 문제에서는 어떤 큰 조합론적 값을 소수 모듈러(대개 $10^9 + 7$)로 계산하도록 요구하는 경우가 많습니다. 이렇게 하는 이유는 참가자들이 큰 수 연산을 사용하지 않고 문제 자체에 집중할 수 있도록 하기 위해서입니다.

일반적인 산술 연산은 모듈러에서도 그리 어렵지 않게 수행할 수 있습니다. 모듈러를 취하고 오버플로에만 주의하면 됩니다. 예를 들면 다음과 같습니다.

c = (a + b) % mod;
c = (a - b + mod) % mod;
c = a * b % mod;

하지만 나눗셈에서는 문제가 생깁니다. 그냥 나누기만 해서는 안 됩니다.

예를 들어 $\frac{8}{2} = 4$이지만

$$ \frac{8 \bmod 5}{2 \bmod 5} = \frac{3}{2} \neq 4 $$

$\frac{1}{a} = a^{-1}$처럼 동작하는 어떤 원소를 찾아서, «나눗셈» 대신 그 원소를 곱해야 합니다. 이러한 원소를 모듈로 $m$에 대한 역원이라고 합니다. 일반적인 나눗셈과 마찬가지로, $a = 0$에 대해서는 모듈로 역원이 정의되지 않습니다.

#이진 거듭제곱을 이용한 방법

페르마의 소정리에 따르면, 임의의 소수 $p$와 임의의 정수 $a$에 대하여 다음이 성립합니다.

$$ a^p \equiv a \pmod p $$ 이제 이 잘 알려진 결과를 $a$로 두 번 «나눕시다»: $$ a^p \equiv a \implies a^{p-1} \equiv 1 \implies a^{p-2} \equiv a^{-1} $$

결과적으로 $a^{p-2}$는 모듈러 곱셈에 대해 $a^{-1}$처럼 동작하며, 이것이 바로 우리에게 필요한 것입니다.

$a^{p-2}$는 이진 거듭제곱으로 $O(\log p)$ 시간에 계산할 수 있습니다.

const int mod = 1e9 + 7;

// бинарное возведение в степень по модулю
int binpow(int a, int n) {
    int res = 1;
    while (n != 0) {
        if (n & 1)
            res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

// находит обратный элемент как a^(p-2)
int inv(int x) {
    return binpow(x, mod - 2);
}

이 방법은 간단하고 빠르지만, 소수 모듈러에 대해서만 동작한다는 점을 기억해야 합니다.

합성수 모듈러의 경우, 오일러 정리에 따라 수 $a$를 $(\phi(m)-1)$ 제곱해야 하며, 이를 위해서는 소인수분해를 구해야 합니다.

#확장 유클리드 알고리즘을 이용하는 방법

확장 유클리드 알고리즘을 사용하면 다음과 같은 형태의 방정식을 정수 범위에서 풀 수 있습니다.

$$ Ax + By = 1 $$ $A$와 $B$ 자리에 각각 $a$와 $m$을 대입합시다: $$ ax + my = 1 $$ 방정식의 해 중 하나가 바로 $a^{-1}$입니다. 왜냐하면 방정식을 $m$으로 나눈 나머지를 취하면 다음을 얻기 때문입니다. $$ ax + my = 1 \iff ax \equiv 1 \iff x \equiv a^{-1} \pmod m $$

거듭제곱을 이용하는 방법에 비해 이 방법이 가지는 장점은 다음과 같습니다.

역원이 존재한다면, 모듈러가 소수가 아니더라도 역원을 찾을 수 있습니다.
이 알고리즘은 손으로 수행하기가 더 쉽습니다.
최적화하면 알고리즘이 조금 더 빠릅니다

하지만 저자 개인적으로는 거의 항상 거듭제곱을 사용합니다.

#간소화된 구현

먼저 구현을 제시하고, 그다음에 이것이 왜 동작하는지 이해해 봅시다:

int inv(int a, int m) {
    if (a == 1)
        return 1;
    return (1 - 1ll * inv(m % a, a) * m) / a + m;
}

이 함수가 실제로 모듈러 역원을 반환한다는 것을 귀납법으로 증명합시다.

기저 사례는 자명합니다: $1 \cdot 1 \equiv 1$.

두 번째 경우에는 공식이 올바른지 확인해 봅시다:

$f(m \bmod a, a) \equiv m^{-1} \pmod a$이므로 $(1 - f(m \bmod a, a) \cdot m)$은 $a$로 나누어떨어집니다.
$\frac{f(m \bmod a, a) \cdot m}{a}$는 $m$으로 나누어떨어지므로, 최종 식은 $m$을 법으로 $\frac{1}{a} = a^{-1}$과 합동입니다.

답이 $0$부터 $(m - 1)$까지의 범위에 들어오는 이유는 독자의 연습 문제로 남겨 두겠습니다.

#모듈러 역원 전처리

대부분의 경우 모듈러 역원은 조합론의 맥락에서 구해야 합니다.

예를 들어, 특히 이항 계수를 계산해야 하는 경우가 많은데, 이를 위해서는 다시 팩토리얼의 역원을 구할 수 있어야 합니다:

$$ C_n^k = \frac{n!}{(n-k)! k!} $$

간단한 방법은 일반 팩토리얼을 미리 계산해 두고 매번 inv를 한두 번 호출하는 것입니다:

int t[maxn]; // факториалы, можно предподсчитать простым циклом

int c(int n, int k) {
    return t[n] * inv(t[k]) % mod * inv(t[n - k]) % mod;
}

// или, почти в два раза быстрее:
int c(int n, int k) {
    return t[n] * inv(t[k] * t[n - k] % mod) % mod;
}

그러나 이 방법은 어떤 조합 공식이 핫 루프 안에 놓이는, 드물지 않은 경우에 점근적 복잡도에 불필요한 로그를 하나 더 추가하게 됩니다. 따라서 자주 쓰이는 모듈러 역원도 미리 계산해 두는 것이 좋습니다.

#역팩토리얼

이미 inv를 작성해 두었다면, $(a!)^{-1}$을 계산하느라 $O(\log m)$번의 연산을 추가로 쓰는 것은 아깝지 않습니다.

그 후 $(a-1)!$의 역원은 다음 식으로 $O(1)$에 계산할 수 있습니다:

$$ (a-1)!^{-1} = (a!)^{-1} \cdot a \equiv \frac{1}{1 \cdot 2 \cdot \ldots \cdot (a-1)} \pmod p $$

나머지 모든 역팩토리얼도 같은 방식으로 이전 값으로부터 반복적으로 계산할 수 있습니다.

// обычные факториалы:
int f[maxn];

f[0] = 1;
for (int i = 1; i < maxn; i++)
    f[i] = i * f[i - 1] % mod;

// обратные факториалы:
int r[maxn];

r[maxn - 1] = inv(f[maxn - 1])
for (int i = maxn - 1; i >= 1; i--)
    r[i-1] = r[i] * i % mod;

또한 $1$부터 $(p - 1)$까지의 모든 수에 대한 역원을 구하는 방법도 있지만, 보통 모듈러가 크기 때문에 자주 쓰이지는 않습니다.

#왜 $10^9+7$입니까?

몇 가지 이유가 있습니다:

이 식은 입력하기가 꽤 쉽습니다(1e9+7).
소수입니다.
충분히 큰 수입니다.
int는 덧셈할 때 오버플로가 발생하지 않습니다.
long long은 곱셈 시 오버플로가 발생하지 않습니다.

참고로, $10^9 + 9$도 이와 똑같은 성질을 모두 가지고 있습니다. 때로는 이것을 사용하기도 합니다.

때로는 $998244353$도 볼 수 있습니다. 이 수는 첫 번째를 제외한 모든 성질을 갖추고 있으며, 그 대신 고속 푸리에 변환의 변형 중 하나에 활용됩니다. 참가자들에게 문제의 주제를 드러내지 않으려고, 이 수와 관련이 없는 문제에도 일부러 이 수를 쓰는 경우가 있습니다.