BOJ 11401 - 이항 계수 3

문제

• \(\dbinom{N}{K}, 1\le N \le 4,000,000, \; 0 \le K \le N\)을 1000000007로 나눈 나머지를 구하는 문제이다.

풀이

• MOD가 \(10^9-7\)이기 때문에, long long으로 충분하다.
• 계산해야 하는 식은 \(\frac{N!}{K!(N-K)!}\)인데, 모듈러 연산은 나눗셈에 대해 분배법칙이 성립하지 않는다.
  따라서 나누려는 숫자의 모듈러 연산에 대한 곱셈의 역원(MMI, Modular Multiplicative Inverse)을 찾고, 그것을 곱해야 한다.
• 이때 확장 유클리드 알고리즘(EEA) 또는 페르마 소정리(FlT)를 사용할 수 있다.
  각 방법을 사용하는 풀이를 순서대로 설명하겠다.
• 그 이전에, 위 조합의 식을 아래와 같이 정리할 수 있다.

\[\begin{aligned} \dbinom{N}{K}=\frac{N!}{K!(N-K)!}&=\frac{K+1\times\dots\times{N}}{1\times\dots\times{N-K}}&&=\prod_{i=0}^{N-K-1}\frac{k+i+1}{i+1}\\ &=\frac{N-K+1\times\dots\times{N}}{1\times\dots\times{K}}&&=\prod_{i=0}^{k-1}\frac{N-i}{K-i} \end{aligned}\]

• 아래 식이 더 간단하므로, 아래 코드에서는 결과의 분자를 a, 분모를 b로 계산하겠다.

EEA

• \(10^9-7\)을 \(M\)이라 하면, 우리는 \(b\cdot b^{-1} \equiv 1 \bmod{M}\)을 만족하는 \(b^{-1}\)의 값을 찾아야 한다.

EEA를 사용하면, 임의의 정수 \(a, b\)에 대해 \(sa+tb= \text{gcd(a, b)}\)를 만족하는 \(s, t\)를 구할 수 있다.
\(a, b\)에 각각 b, \(M\)을 넣으면, \(sb+tM=\text{gcd(b, M)}\)을 만족하는 \(s, t\)를 얻을 수 있다.
이때 \(M\)이 소수이기 때문에, \(gcd(b, M)= 1\)이다.
따라서 \(sb+tM\equiv sb \equiv 1 \bmod{M}\)이다.
즉, \(s\)가 b의 MMI인 것이다.

• EEA의 동작 과정은 아래 코드의 modinv 함수를 읽는 것이 더 직관적일 것이다.
• 초기화는 modinv(a, b, 1, 0, 0, 1)으로 실행해야 한다.

EEA 코드

#include <iostream>
#include <tuple>
#include <vector>
#define v vector
#define MOD 1000000007

using namespace std;
using vi = vector<int>;
using ll = long long;

tuple<ll, ll> modinv(ll a, ll b, ll s1, ll s2, ll t1, ll t2) {
    if (b == 0) return {s1, t1};
    ll q = a / b, t;
    t = s2;
    s2 = s1 - q * s2;
    s1 = t;
    t = t2;
    t2 = t1 - q * t2;
    t1 = t;
    return modinv(b, a % b, s1, s2, t1, t2);
}

ll getinv(ll a, ll b) {
    ll f = 0;
    if (a < b) {
        ll t = a;
        a = b;
        b = t;
        f = 1;
    }
    if (f)
        f = get<1>(modinv(a, b, 1, 0, 0, 1));
    else
        f = get<0>(modinv(a, b, 1, 0, 0, 1));
    if (f < 0) f += MOD;
    return f;
}

int main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    ll n, k, a = 1, b = 1, ans;
    cin >> n >> k;
    for (int i = 0; i < k; i++) {
        a = a * (n - i) % MOD;
        b = b * (k - i) % MOD;
    }
    ans = a * getinv(b, MOD) % MOD;
    cout << ans;
    return 0;
}

FlT

\[\begin{aligned} p\in\mathbb{P},\;a\in\mathbb{Z}&\implies a^p\equiv a\bmod{p}\\ p\in\mathbb{P},\;p\perp a\in\mathbb{Z}&\implies a^{p-1}\equiv 1\bmod{p} \end{aligned}\]

• 이 문제에서는 \(M\)이 소수이기 때문에, \(b^{M-1}\equiv b\cdot b^{M-2} \equiv 1\bmod M\)이므로, \(b^{M-2}\)가 MMI이다.
• 하지만 \(M\)이 매우 크기 때문에, 아래 코드에서의 pow 함수와 같이 Exponentiation By Squaring을 사용해야 한다.

FlT 코드

#include <iostream>
#include <vector>
#define v vector
#define MOD 1000000007

using namespace std;
using vi = vector<int>;
using ll = long long;


ll pow(ll x, ll p) {
    ll r = 1;
    while (p) {
        if (p & 0x01) r = r * x % MOD;
        x = x * x % MOD;
        p >>= 1;
    }
    return r;
}

int main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    ll n, k, a = 1, b = 1, ans;
    cin >> n >> k;
    for (int i = 0; i < k; i++) {
        a = a * (n - i) % MOD;
        b = b * (k - i) % MOD;
    }
    ans = a * pow(b, MOD - 2) % MOD;
    cout << ans;
    return 0;
}

풀고나서

• 아래 순서대로 코드를 개선했다.
  1. 팩토리얼 배열 직접 계산 + EEA
  2. 팩토리얼 배열 직접 계산 + FlT
  3. 조합 항 정리(O(k)) + FlT
• 1 -> 2는 공간복잡도가 줄어들었고, 2 -> 3은 시간복잡도가 줄어들었다.
• ps에서 자주 사용되는 정수론 지식을 정리할 수 있어서 좋았다.
  모듈러 연산부터 FlT까지의 내용은 따로 포스팅해야겠다.