PS를 위한 정수론 - (4) 이항 계수 (nCr mod P) 구하는 다양한 방법

 

 

이항 계수 $_{n}C_r$을 소수 $p$로 나눈 나머지를 빠르게 구하는 다양한 방법들을 알아보자.

 

기본적으로 특정 $n$과 $r$에 대해서 이항 계수 $_{n}C_r$을 구하는 시간은 $O(1)$이 되어야 할 때,

즉, 많은 쿼리가 들어와도 문제가 없는 경우를 고려한다. 

지금부터 소개하는 방법들의 시간복잡도는 전처리하는데 필요한 시간을 의미한다. 

 

너무 작은 $n$과 $r$에 대해서 직접 분자와 분모를 모두 곱하는 경우는 생략한다. 

 

  1. $O(n^2)$

 

www.acmicpc.net/problem/11051

11051번: 이항 계수 2

 

파스칼의 삼각형을 이용하는 방식이다. 

이항 계수에 대한 성질로 매우 잘 알려져있는 $_{n}C_r = \;_{n-1}C_{r-1} +\; _{n-1}C_r$ 을 이용해서 dp table을 만든다. 

따라서 이를 전처리하기 위해서는 $O(n^2)$에 해당하는 시간과 메모리가 필요하다. 

 

가장 단순한 방법이지만, $n$과 $r$이 10000 이상 커진다면 이용하기 힘든 방식이다.

물론 $n$과 $r$이 충분히 작다면 이 방식을 이용하지 않을 이유가 없다. 

 

[소스 코드]

#include<bits/stdc++.h>
using namespace std;
using ll = long long;
const int MOD = (int)1e9 + 7;

ll C[1000][1000];
int main(void) {
	int n, r;
	for (int i = 1; i <= 1000; i++) {
		for (int j = 0; j <= i; j++) {
			if (j == 0 || j == i) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % MOD;
		}
	}
	cout << C[n][r];
}

 

 

 

  2. $O(nlogp)$

 

www.acmicpc.net/problem/11401

11401번: 이항 계수 3

이처럼 $n$과 $r$이 400만인 경우는 어떻게 구할까?

위의 방식으로는 400만 x 400만의 배열 크기와 시간이 필요하므로 불가능하다는 것을 당연히 알 수 있다. 

 

페르마의 소정리를 떠올려보자. (링크 : rebro.kr/105)

n이 10000보다 커진다면, 파스칼의 삼각형을 이용하는 방식을 사용할 수 없으므로

$_{n}C_r = n! / (r!(n-r)!) $식을 이용해서 분모와 분자를 직접 계산하는 방식을 이용해줘야 한다. 

다만, 우리가 구하는 이항 계수는 $p$에 대한 나머지를 구하는데, 분모가 존재하므로 계산하기가 쉽지 않다. 

따라서 페르마의 소정리인 $a^{p-1} \equiv 1\;(mod\;p)$를 이용한다. 

이 식을 이용하면 $a^{p-2} \equiv a^{-1}\;(mod\;p)$, 즉 $a$의 $mod\;p$에 대한 곱셈의 역원이 $a^{p-2}$인 것을 알 수 있다. 

 

따라서 $mod\;p$에 대해서 $r!(n-r)!$을 나누는 것이 아니라, $r!(n-r)!$의 역원인 $(r!(n-r)!)^{p-2}$을 곱하는 방식이다. 

 

$a^k$를 구하는데 필요한 시간은 $O(logk)$이면 충분하다. "분할 정복을 이용한 거듭제곱" 방식이다. 

혹시 이를 모른다면 꼭 공부를 하는 것을 권장한다. 

 

예를 들어서 $a^{58}$을 계산해야 한다고 가정하자. $58$을 이진수로 나타내면 $111010$로 나타난다.

따라서 $a^{58}$의 지수를 $2^5 + 2^4 + 2^3 + 2^1$로 표현하면 $a^{58} = a^{32}\times a^{16}\times a^8 \times a^2$와 같이 표현이 가능하고, 

이진수의 오른쪽 끝부터 시작해서 매번 a를 제곱시켜가면서, $1$인 자릿수인 경우에는 결괏값에 곱해주면 된다. 

 

여기까지는 하나의 $_{n}C_r$을 구하는 과정이고, 최종적으로 우리가 원하는 경우는 많은 쿼리가 들어와도 문제가 없는 경우이다. 

www.acmicpc.net/problem/13977

13977번: 이항 계수와 쿼리

매번 $_{n}C_r$을 $O(1)$만에 구하기 위해서는 팩토리얼(factorial) 값에 대해서 미리 전처리를 해두어야 한다. 

$1!$부터 $n!$까지 배열에 계산을 해두고, 각 팩토리얼값마다 역원을 계산$((O(logp))$해두면 $_{n}C_r$을 $O(1)$만에 계산할 수 있다. 

 

코드는 3번에서 확인하자. 

 

  3. $O(n + logp)$

 

2번 방식에서 약간만 변형시킨 것이다. 

 

원래는 각 팩토리얼값마다 역원을 각각 계산해주었는데, 잘 생각해보면 팩토리얼의 역원에도 점화식이 생긴다. 

$((n-1)!)^{-1} = n*(n!)^{-1}$ 을 만족하기 때문에, 처음에 $1!$부터 $n!$까지의 값과 $n!$의 역원을 미리 계산해두면

$1!$부터 $(n-1)!$의 역원까지는 모두 $O(1)$에 계산할 수 있다. 

 

사실 작정하고 이항계수에 대해서 낸 문제가 아니라면 기본적으로 N이 1000만 이상 주어지지 않기 때문에 보통 대부분의 문제들은 여기까지만 알아도 다 해결할 수 있다.

 

 

[소스 코드]

#include<bits/stdc++.h>
using namespace std;
using ll = long long;
const int MOD = (int)1e9 + 7;

#define MN 4000000
ll fac[MN+10], facinv[MN+10];
ll mpow(ll a, ll x) {
	ll res = 1;
	while (x > 0) {
		if (x % 2) res = (res*a) % MOD;
		a = (a*a) % MOD;
		x /= 2;
	}
	return res;
}
void fac_init() {
	fac[0] = 1;
	for (int i = 1; i <= MN; i++) fac[i] = fac[i - 1] * i % MOD;
	facinv[MN] = mpow(fac[MN], MOD - 2);
	for (int i = MN - 1; i >= 0; i--) facinv[i] = facinv[i + 1] * (i + 1) % MOD;
}
ll C(ll n, ll r) {
	return ((fac[n] * facinv[r]) % MOD) * facinv[n - r] % MOD;
}
int main(void) {
	ll n, r;
	fac_init();
	cout << C(n, r);
}

 

 

 

  4. $O(n)$

 

$O(n + logp)$와 비교해서 시간적인 부분은 거의 동일하다. 

(실제로 코드 구현에 있어서는 추가적인 연산들로 인해 $O(n + logp)$의 시간이 덜 소모된다) 

다만 구하는 과정이 짧고 꽤나 참신해서 살펴볼만하다. 

 

3번에서 $O(logp)$가 필요한 이유는 처음에 $n!$에 대한 역원을 직접 구해줘야 했기 때문이다. 

이를 어떻게 $O(1)$만에 구할 수 있을까?

 

정수 $k$에 대한 역원은 DP로 간단하게 구할 수 있고, 점화식은 다음과 같다. 

$1^{-1} \equiv 1$이고, $k^{-1} \equiv -(p / k) \times (p\%k)^{-1} \; (mod\;p)$ 

 

증명 과정을 살펴보자.

 

먼저 $p = (p/k)*k + p\%k$ 를 만족한다.

양변에 $mod\;p$를 씌워주면 $0 \equiv (p/k)*k + p\%k\;(mod\;p)$  ,  $(-p/k)*k \equiv p\%k \;(mod\;p)$ 가 된다. 

 

따라서 $((-p/k)*k)^{-1} \equiv (p\%k)^{-1} \;(mod\;p)$ 를 만족하고,

양변에 $(-p/k)$를 곱해주면 최종적으로 $k^{-1} \equiv (p\%k)^{-1} \times (-p/k)$ 가 나온다. 

 

여기서 나오는 $(p/k)$는 정수이고, $p\%k$는 항상 $k$보다 작기 때문에 이전의 dp값으로 계산해줄 수 있다. 

 

[소스 코드]

#include<bits/stdc++.h>
using namespace std;
using ll = long long;
const int MOD = (int)1e9 + 7;

#define MN 4000000
ll fac[MN+10], facinv[MN+10], inv[MN+10];
void fac_init() {
	fac[0] = facinv[0] = inv[1] = 1;
	for (int i = 1; i <= MN; i++) {
		inv[i+1] = (MOD - MOD / (i+1)) * inv[MOD%(i+1)] % MOD;
		fac[i] = fac[i - 1] * i % MOD;
		facinv[i] = facinv[i-1]* inv[i] % MOD;
	}
}
ll C(ll n, ll r) {
	return ((fac[n] * facinv[r]) % MOD) * facinv[n - r] % MOD;
}
int main(void) {
	ios::sync_with_stdio(false); cin.tie(nullptr); cout.tie(nullptr);
	fac_init();
	ll n, k; cin >> n >> k;
    cout << C(n, k) ;
}

 

 

  5. $O(p)$

 

www.acmicpc.net/problem/11402

11402번: 이항 계수 4

이처럼 $n$이 매우 큰 경우에는 위에서 설명한 어떤 방법으로도 $O(n)$보다 작게 계산할 수 없다. 

이때, Lucas Theorem (뤼카의 정리)이라는 정리를 사용할 수 있다. 

 

뤼카의 정리는 이처럼 $n$이 크고 $p$가 작은 경우에 이항계수를 $O(p)$만에 구할 수 있는 방법이다. 

 

- Lucas Theorem (뤼카의 정리)

 

음이 아닌 정수 $n$과 $r$, 소수 $p$에 대해서 $n$과 $r$를 $p$진법으로 나타내면 다음과 같다.

 

$n = n_kp^k + n_{k-1}p^{k-1} + n_{k-2}p^{k-2} + ... + n_1p + n_0$

$r = r_kp^k + r_{k-1}p^{k-1} + r_{k-2}p^{k-2} + ... + r_1p + r_0$

 

이때, $_{n}C_r \equiv \prod_{i=0}^{k} {_{n_i}C_{r_i}} \;\; (mod\;p)$ 가 성립한다는 것이 뤼카의 정리이다. 

---

증명과정은 다음과 같다. 

 

$\sum_{r=0}^{n} {_{n}C_r} x^r $ 의 전개과정을 살펴보자. 

$\sum_{r=0}^{n} {_{n}C_r} x^r \equiv (1+x)^n \equiv (1+x)^{n_k p^k + n_{k-1}p^{k-1} + ... + n_1p + n_0} \equiv \prod_{i=0}^{k} [(1+x)^{p^i}]^{n_i}$를 만족한다.

 

이때, $(1+x)^{p^n}$은 이항정리에 의해서 $_{p^n}C_0 x^0 + _{p^n}C_1 x^1 + ... + _{p^n}C_{p^n} x^{p^n}$ 을 만족하고, $p$가 소수이므로 모든 $1\leq i \leq p^n-1$에 대해서 $_{p^n}C_i$는 $p$로 나누어 떨어진다.

따라서, $(1+x)^{p^n} \equiv 1+x^{p^n} \;\;(mod\;p)$를 만족한다. 

 

그러므로 $\prod_{i=0}^{k} [(1+x)^{p^i}]^{n_i} \equiv \prod_{i=0}^{k} [1+x^{p^i}]^{n_i} \;\; (mod \;p)$ 로 만들어진다. 

 

여기서 다시 안의 식에 이항 정리를 이용하면 $\prod_{i=0}^{k} \,[\,\sum_{r_i = 0}^{n_i} \;{_{n_i}C_{r_i} x^{r_ip^i}}\,]$가 되고,

$[ \,\,]$ 안의 식을 전개한 후 정리하면 $\sum_{r=0}^{n} [\,(\prod_{i=0}^{k} {_{n_i}C_{r_i}}) x^r\,]$ 로 만들어진다. 

 

따라서, $_{n}C_r \equiv \prod_{i=0}^{k} {_{n_i}C_{r_i}} \;\; (mod \;p)$ 가 성립한다. 

 

증명 과정은 복잡하지만, 정리를 코드로 구현하면 생각보다 간단하다. 

쿼리당 걸리는 시간은 $O(logn \;/\; logp)$ 이므로 상수와 다름없게 생각해도 된다. 

 

[소스 코드]

#include<bits/stdc++.h>
#define MAX_P 2001
using namespace std;
using ll = long long;
const int MOD = 7;
ll fac[MAX_P];

ll mpow(ll a, ll k, ll mod) {
	ll res = 1;
	while (k > 0) {
		if (k % 2) res = (res * a) % mod;
		a = (a*a) % mod;
		k /= 2;
	}
	return res;
}

int main(void) {
	ll n, r; 
	ll ans = 1;
	while (n || r) {
		int N = n % MOD;
		int R = r % MOD;
        if(R > N) {
            ans = 0;
            break;
        }
		ans *= fac[N] * mpow(fac[R] * fac[N - R], MOD - 2, MOD) % MOD;
		ans %= MOD;
		n /= MOD;
		r /= MOD;
	}
	cout << ans;
}

 

 

 

  6. 연습 문제

 

[BOJ 13977] 이항 계수와 쿼리  Gold I

 

[BOJ 11402] 이항 계수 4  Platinum V

 

[BOJ 15791] 세진이의 미팅  Gold I

 

[BOJ 16134] 조합  Gold I

 

[BOJ 20296] 폰친구  Platinum IV

 

[Codeforces 1462E2] Close Tuples (Hard version)

 

[Codeforces 1445D] Divide and Sum

 

 

 

본 글은 구사과님의 블로그 (koosaga.com/63)에서 영감을 받은 글입니다. 

 

PC로 보시는 것을 권장합니다. 

피드백은 언제나 환영입니다. 댓글로 달아주세요 ^-^