이항 계수(Binomial Coefficient)를 구하는 다양한 알고리즘

이항 계수를 효율적으로 구하는 알고리즘을 정리한 글입니다.
KKH
·
2023년 04월 05일

들어가며

조합론에서, 이항 계수(Binomial Coefficient)란 이항식을 이항 정리로 전개했을 때 각 항의 계수이며, 주어진 크기의 (순서 없는) 조합의 가짓수라고 말합니다.

이렇게 들으면 이해가 안 되실 수 있으니, 학교에서 배운 공식과 기호를 통해 나타내 보겠습니다.

정의

자연수 및 정수 nn과 kk가 주어졌을 때, 이항 계수 (nk)\binom{n}{k}는 아래와 같습니다.


(nk)={n!k!⋅(n−k)!(0≤k≤n)0(k<0)0(k>n)\binom{n}{k} = \left\{ \begin{array}{ll} \frac{n!}{k! \cdot (n - k)!} & (0 \leq k \leq n) \\ 0 & (k < 0) \\ 0 & (k > n) \end{array} \right.

또한, 이항 계수 (nk)\binom{n}{k}는 nCknCk이나 C(n,k)C(n, k)로 나타내기도 합니다.

이 글에서는, 이항 계수를 구하는 기초적인 알고리즘부터, 매우 큰 이항 계수를 빠르게 구하는 알고리즘까지 정리해 보았습니다.

1. Divide and Conquer

이항 계수는 아래와 같은 두 가지의 성질이 있습니다.


(nk)=(n−1k)+(n−1k−1)\binom{n}{k} = \binom{n - 1}{k} + \binom{n - 1}{k - 1}

(n0)=(nn)=1\binom{n}{0} = \binom{n}{n} = 1

따라서 저희는 위와 같은 두 가지의 성질을 이용해 재귀적으로 분할 정복을 수행하는 bino_dnc 함수를 구현할 수 있습니다. 예시 코드는 아래와 같습니다.

#include <bits/stdc++.h>

using namespace std;

int bino_dnc(int n, int k) {
    if (k == 0 || k == n) return 1;
    return bino_dnc(n - 1, k) + bino_dnc(n - 1, k - 1);
}

int main() {
    int n, k;
    cin >> n >> k;
    cout << bino_dnc(n, k);

    return 0;
}

분할 정복으로 이항 계수를 구할 때의 시간 복잡도는 O(2n),O(2n/n)O(2^n), O(2^n / \sqrt{n})로 근사할 수 있습니다. (증명은 어려우니 생략... - 참고 링크)

구현하기에는 되게 쉽고 직관적이지만, 위의 알고리즘에는 치명적인 단점이 있습니다. 바로 시간이 너무 오래 걸리고, 중복된 계산이 너무 많다는 점인데요. 이를 개선하는 알고리즘을 알아봅시다.

2. Dynamic Programming

Dynamic Programming(동적 계획법)은 큰 문제를 부분 문제로 쪼개어 부분 문제를 해결한 뒤, 그 결과를 통해 큰 문제를 해결해 나가는 알고리즘 기법입니다.

큰 문제를 쪼개어 해결하는 것에 있어서는 분할 정복과 비슷해 보일 수 있지만, 분할 정복과 동적 계획법의 가장 큰 차이는 Memoization에 있습니다. Memoization을 통한 최적화로 불필요한 중복 계산을 줄일 수 있습니다.

앞서 설명한 이항 계수의 성질을 활용해 아래와 같은 점화식을 세울 수 있습니다.


DP[i][j]={DP[i−1][j]+DP[i−1][j−1](if 0<j<i)1(if j=0 or j=i)DP[i][j] = \left\{ \begin{array}{ll} DP[i - 1][j] + DP[i - 1][j - 1] & (\text{if } 0 < j < i) \\ 1 & (\text{if } j=0 \text{ or } j=i) \end{array} \right.

그리고 위의 점화식을 토대로 bino_dp 함수를 구현한 예시 코드는 아래와 같습니다.

#include <bits/stdc++.h>

using namespace std;

const int MOD = 1e9 + 7;
int DP[1001][1001];

int bino_dp(int n, int k) {
    for (int i = 1; i <= 1000; i++) {
        for (int j = 0; j <= min(i, k); j++) {
            if (j == 0 || j == i) DP[i][j] = 1;
            else DP[i][j] = (DP[i - 1][j] + DP[i - 1][j - 1]) % MOD;
        }
    }
    return DP[n][k];
}

int main() {
    int n, k;
    cin >> n >> k;
    cout << bino_dp(n, k);

    return 0;
}

알고리즘 문제를 해결할 때 정답이 자료형의 표현 가능 범위를 벗어나게 되는 경우, 정답에 대해 적당히 큰 소수의 모듈러 연산을 수행한 값을 구해야 할 때가 많습니다.

앞으로 설명드릴 코드에서는 이항 계수의 NN과 KK의 값이 커질 때, 결과 값이 컴퓨터가 표현 가능한 범위에서 벗어나게 되므로 모듈러 연산한 값을 결과로 반환합니다. 이처럼 특정 문제의 결과 값이 매우 커질 때, 모듈러 연산이 자주 활용되기 때문에 이에 대한 성질은 알고 계시는 것이 좋습니다.

동적 계획법의 시간 복잡도를 계산해 보면, O(NK)O(NK)가 되는 것을 확인하실 수 있을 겁니다. 이는 분할 정복보다 확실히 더 빠른 시간 내에 수행할 수 있습니다.

하지만, 이 알고리즘에도 치명적인 단점이 있습니다... 바로 NN과 KK의 값이 10,00010,000 이상 커진다면 컴퓨터가 수행할 수 없을 것입니다. NN과 KK의 제한이 작다면, 고려해볼 만한 방법 중 하나입니다. 그럼, 이보다 더 나은 방법을 알아보도록 하겠습니다.

3. 페르마의 소정리

이항 계수에서 NN과 KK의 값이 커진다면, 앞서 설명드린 방법으로 구하기가 거의 불가능할 것입니다. 이 때, 페르마의 소정리를 통해 빠르게 이항 계수를 구할 수 있습니다.

우선, 페르마의 소정리에 대해 알아보겠습니다. 페르마의 소정리란,

  1. pp가 소수이면, 모든 정수 aa에 대해
ap≡a(modp)a^p \equiv a \pmod{p}
  1. pp가 소수이고 aa가 pp의 배수가 아니면
ap−1≡1(modp)a^{p-1} \equiv 1 \pmod{p}

위 두 정리는 동치인 정리입니다.

우리는 이항 계수를 구할 때 위의 정리를 활용할 수 있습니다. 우선, 이항 계수의 공식을 다시 한 번 봅시다.

(nk)=n!k!⋅(n−k)!\binom{n}{k} = \frac{n!}{k! \cdot (n - k)!}

우리는 위 식을 이용하여 분모와 분자를 직접 계산하는 방식으로 구해야 합니다. 팩토리얼은 그냥 계산하면 되고, 문제는 분모입니다. 분모 때문에 모듈러 연산을 하기가 쉽지 않습니다. 이 때, 페르마의 소정리를 이용합니다.

페르마의 소정리에서, 아래와 같은 식을 도출해낼 수 있습니다.

ap−1≡1(modp)a^{p-1} \equiv 1 \pmod{p}
⇒a⋅ap−2≡1(modp)\Rightarrow a \cdot a^{p-2} \equiv 1 \pmod{p}
⇒ap−2≡a−1(modp)\Rightarrow a^{p-2} \equiv a^{-1} \pmod{p}

즉, 우리는 페르마의 소정리를 통해 간단하게 분모의 곱셈 역원을 구할 수 있고, 이를 곱해준다면 분모를 나눠주는 것과 같은 결과를 얻을 수 있습니다.

정리해 보면,

(nk)(modp)\binom{n}{k} \pmod{p}
=n!k!⋅(n−k)!(modp)= \frac{n!}{k! \cdot (n - k)!} \pmod{p}
=n!⋅(k!⋅(n−k)!)p−2(modp)= n! \cdot (k! \cdot (n - k)!)^{p-2} \pmod{p}
=n!⋅(k!)p−2⋅((n−k)!)p−2(modp)= n! \cdot (k!)^{p-2} \cdot ((n - k)!)^{p-2} \pmod{p}

가 됩니다. 우리는 분할 정복을 이용한 거듭제곱 알고리즘을 활용해 ana^n을 O(logn)O(logn) 안에 구할 수 있습니다. 팩토리얼 값은 전처리를 통해 1!∼N!1! \sim N!까지 미리 구해놓는다면, 우리는 팩토리얼의 역원을 구해 곱하기만 하면 되니 알고리즘의 시간복잡도는 O(n)O(n)이 되는 것을 알 수 있습니다.

아래는 페르마의 소정리를 활용한 이항 계수를 구하는 예시 코드입니다.

#include <bits/stdc++.h>
#define MAX_N 4000000

using namespace std;
using ll = long long;

const int MOD = 1e9 + 7;
ll fac[MAX_N + 1];

ll pow_dnc(int a, int n) {
    if (n == 0) return 1 % MOD;
    ll u = pow_dnc(a, n / 2);
    u = (u * u) % MOD;
    if (n % 2 == 1) u = (u * a) % MOD;
    return u;
}

void fac_init() {
	fac[0] = 1;
	for (int i = 1; i <= MAX_N; i++) fac[i] = fac[i - 1] * i % MOD;
}

ll bino_fermat(int n, int k) {
	return ((fac[n] * pow_dnc(fac[k], MOD - 2)) % MOD) * pow_dnc(fac[n - k], MOD - 2) % MOD;
}

int main() {
    int n, k;
    cin >> n >> k;
    fac_init();
    cout << bino_fermat(n, k);

    return 0;
}

위 코드의 시간 복잡도는 O(n)+O(logn)=O(n)O(n) + O(logn) = O(n)이 됩니다. 이정도면 충분히 큰 NN에 대해서도 빠르게 구할 수 있습니다. 페르마의 소정리는 다양한 곳에서 활용되기 때문에 알아두시면 좋습니다.

최적화

위에서 소개드린 코드에서는, 각 팩토리얼 값마다 역원을 계산해 주는 방식이었습니다. 하지만, 팩토리얼의 식을 잘 보면 팩토리얼의 역원에도 점화식이 생기는 것을 볼 수 있습니다.

((n−1)!)−1=n⋅1n!=n⋅(n!)−1((n - 1)!)^{-1} = n \cdot \frac{1}{n!} = n \cdot (n!)^{-1}

이 성립하기 때문에, 전처리 단계에서 1!∼N!1! \sim N!의 값과 N!N!의 역원을 구해준 뒤, (N−1)!∼1!(N - 1)! \sim 1!까지의 역원을 구해주면 이항 계수의 계산을 O(1)O(1)에 처리할 수 있습니다. 예시 코드는 아래와 같습니다.

#include <bits/stdc++.h>
#define MAX_N 4000000

using namespace std;
using ll = long long;

const int MOD = 1e9 + 7;
ll fac[MAX_N + 1], fac_inv[MAX_N + 1];

ll pow_dnc(int a, int n) {
    if (n == 0) return 1 % MOD;
    ll u = pow_dnc(a, n / 2);
    u = (u * u) % MOD;
    if (n % 2 == 1) u = (u * a) % MOD;
    return u;
}

void fac_init() {
	fac[0] = 1;
	for (int i = 1; i <= MAX_N; i++) fac[i] = fac[i - 1] * i % MOD;
	fac_inv[MAX_N] = pow_dnc(fac[MAX_N], MOD - 2);
	for (int i = MAX_N - 1; i >= 0; i--) fac_inv[i] = fac_inv[i + 1] * (i + 1) % MOD;
}

ll bino_fermat_optimize(int n, int k) {
	return ((fac[n] * fac_inv[k]) % MOD) * fac_inv[n - k] % MOD;
}

int main() {
    int n, k;
    cin >> n >> k;
    fac_init();
    cout << bino_fermat_optimize(n, k);

    return 0;
}

4. 뤼카의 정리

만약 NN이 매우 크다면, O(N)O(N)의 시간으로도 부족할 수 있습니다. 이 때에는 더 효율적인 알고리즘을 사용해야 하는데, 그것이 바로 뤼카의 정리입니다.

뤼카의 정리는 NN이 매우 크고, PP가 적당히 작은 수일 때, 이항 계수를 O(P)O(P)에 구할 수 있는 방법입니다.

우선, 음이 아닌 정수 nn과 kk를 pp-진법으로 전개하면 다음과 같습니다.

n=nm⋅pm+nm−1⋅pm−1+nm−2⋅pm−2+⋯+n1⋅p+n0k=km⋅pm+km−1⋅pm−1+km−2⋅pm−2+⋯+k1⋅p+k0n = n_{m} \cdot p^{m} + n_{m-1} \cdot p^{m-1} + n_{m-2} \cdot p^{m-2} + \cdots + n_{1} \cdot p + n_{0} \\ k = k_{m} \cdot p^{m} + k_{m-1} \cdot p^{m-1} + k_{m-2} \cdot p^{m-2} + \cdots + k_{1} \cdot p + k_{0}

이 때, 다음 합동식이 성립합니다.

(nk)≡∏i=0m(niki)(modp)\binom{n}{k} \equiv \prod_{i=0}^{m} \binom{n_{i}}{k_{i}} \pmod{p}

(증명 과정은 저도 아직 이해를 못 했기 때문에... 이해를 하게 된다면 다시 작성하는 걸로...)

아래는 앞서 설명한 뤼카의 정리로 매우 큰 NN과 적당히 작은 PP에 대한 이항 계수를 구하는 예시 코드입니다.

#include <bits/stdc++.h>
#define MAX_P 2000

using namespace std;
using ll = long long;

const int MOD = 7;
ll N, K, fac[MAX_P + 1];
vector<ll> A, B;

ll pow_dnc(int a, int n) {
    if (n == 0) return 1 % MOD;
    ll u = pow_dnc(a, n / 2);
    u = (u * u) % MOD;
    if (n % 2 == 1) u = (u * a) % MOD;
    return u;
}

void fac_init() {
	fac[0] = 1;
	for (int i = 1; i <= MAX_P; i++) fac[i] = fac[i - 1] * i % MOD;
}

void digits_init() {
    while (N > 0) {
        A.push_back(N % MOD);
        N /= MOD;
    }
    while (K > 0) {
        B.push_back(K % MOD);
        K /= MOD;
    }
}

ll bino_fermat(int n, int k) {
	return ((fac[n] * pow_dnc(fac[k], MOD - 2)) % MOD) * pow_dnc(fac[n - k], MOD - 2) % MOD;
}

ll bino_lucas() {
    ll RES = 1;
    for (int i = 0; i < A.size(); i++) {
        int n, k;
        n = A[i];
        if (i >= B.size()) k = 0;
        else k = B[i];
        RES = (RES * bino_fermat(n, k)) % MOD;
    }
	return RES;
}

int main() {
    cin >> N >> K;
	fac_init();
	digits_init();
	cout << bino_lucas();

    return 0;
}

5. 중국인의 나머지 정리 (...?)

이 부분은 저도 이해를 하지 못했습니다. 머리가 더 좋아지면 작성하는 걸로...

결론

이항 계수를 계산하는 많은 알고리즘을 살펴봤습니다. 문제의 조건에 맞춰 적절한 알고리즘을 사용하는 것이 중요할 것입니다.

kbhetrr

Copyright © 2024, All right reserved.