이항 계수 (Binomial Coefficient)
-
$ {\binom {n}{k}}$, ${{n}C_{k}}$, ${C(n,k)}$ 등으로 표기하고, 계산 공식은 0≤k≤n 일 때, $n!/\left(k!(n-k)!\right)$ 이다.
잘 쓰이는 일이 없지만, k>n 일때에는 C(n,k)=0 이다.
k<0일 때는 설명이 다르다
-
C(n,k)의 값은 빠르게 커지므로, n>=35면 32비트 정수형의 범위를 넘어가고 n>=68이면 64비트 정수형의 범위를 넘어가게 된다. 그래서 보통 문제에서는 모듈러를 취한 값을 요구하게 된다.
C(n,k)에서 k>n-k 인 경우에는 C(n,n-k)을 구해도 값이 똑같으므로, 아래부터는 항상 k<n-k임을 가정하고 설명하겠다.
관련 정리
쿠머의 정리
-
어떤 소수 p에 대해서, p^k가 C(n,m)을 나누어 떨어지게 하는 가장 큰 k값을 찾는 문제.
$ \nu _{p}\left({\binom {n}{m}}\right)={\dfrac {S_{p}(m)+S_{p}(n-m)-S_{p}(n)}{p-1}} $
직관적인 활용예로, C(n,m)의 홀짝성을 찾는 문제가 있다.
p=2일때 $\nu_{p}({\binom {n}{m}}) $ 이 0이면 홀수, 1이상이면 짝수라는 뜻이므로, 2진수끼리 더할때 받아올림이 생기는지 없는지만 확인하면 된다. 받아올림이 생기는 경우는 같은 자리에 bit값이 둘다 1인 경우뿐이므로, 결국 이것은 그냥 (n-m)과 m의 bitwise and 연산 결과가 0인지 아닌지만 보면 된다. 0이면 홀수, 1이면 짝수이다.
뤼카의 정리
-
사실상 쿠머의 정리를 조금 다른 형태로 나타낸 정리이다.
$ {\binom {m}{n}}\equiv \prod _{i=0}^{k}{\binom {m_{i}}{n_{i}}}{\pmod {p}} $
이항계수를 소수로 나눈 나머지를 구해야 할때, 더 작은 값에 대한 이항계수들의 곱으로 구할수 있는 방법을 제공한다. 이는 실제로 PS에서 이항계수를 작은 p로 나눈 나머지를 구해야 할때 사용된다.
정확한 값을 계산
직접 구현하는 것도 어렵지는 않지만, 파이썬에서는 이미 구현된 math.comb(n, k) 함수가 있다.
파이썬은 자체적으로 biginteger가 지원되므로 상관 없지만, 다른 일반적인 언어에서 구현할때는 integer overflow에 주의해야 한다.
n!를 먼저 계산해놓고 k!(n-k)! 로 나누려 하면 n!를 계산하는 과정에서 값이 overflow 될 수 있다. n이 20보다 크면, n!이 64비트 정수형의 범위를 초과한다.
C(n,k) = C(n-1,k-1)*n/k 의 형태로 생각해서, 곱하기와 나누기를 번갈아가면서 계산해주면 n⇐62 까지는 64비트 정수형을 갖고 계산 가능하다
ans = 1
for i in range(1, k + 1):
ans = ans * (n - k + i) // i
return ans
파스칼의 삼각형 방식으로 계산하면 속도는 느리지만, 조금 더 큰 n에 대해서까지 계산할 수 있다. 하지만 그래봐야 어차피 68이 넘으면 결과값 자체가 범위 초과이므로 큰 도움은 안되고, 결국 그보다 큰 n에 대해서는 biginteger를 사용해야 한다.
파스칼의 삼각형
이항계수를 큰 소수로 나눈 나머지를 구할 때
값을 한번만 계산할 때
모듈러 나눗셈을 이용해서, 팩토리얼 형태의 공식대로 계산하는게 가장 깔끔하다. C(n,k) = (n-k+1)*(n-k+2)*…*(n-1)*n / (k!) 형태가 되므로 분자, 분모 모두 k번의 곱하기 연산이 필요하고, 분모의 모듈러 역원을 계산하는 데에 O(logP)가 걸리므로 O(k+logP)가 된다.
팩토리얼을
다항식 다중계산을 이용해서 계산한다면 O(k)보다도 빠르게 계산할수도 있긴 하다..
아래에서 언급할 방법을 이용하면, k!의 모듈러 역원을 O(logP)에 계산하지 않고, O(k)에 계산할 수도 있다. 하지만 보통은 logP가 k보다 작아서 무시할수 있고, 게다가 O(k)의 메모리까지 필요로 해서 비효율적.
파이썬에서는 그냥 modular 연산 없이 거대한 이항계수를 math.comb()로 정확히 구하고, 최종적으로 그 값에 modular연산을 딱 한번 적용하는 다소 무식한 방법도 의외로 강력하다. n<5000 정도까지도 위 방법에 비해서 속도가 크게 차이나지 않는다.
코드
관련 문제
값을 여러번 계산해야 할 때
그냥 C(n,k)를 한번만 구하고 끝인 것이 아니라, C(n1,k1), C(n2,k2), …, C(nq,kq) 를 모두 구해야 하는 경우이다 (ni ≤ n)
-
전처리 단계에서 fact[i] = (i! % P) 테이블을 미리 모두 계산해두면, 쿼리 한번에 대해서 C(n,k) mod P = fact[n] * inv(fact[k]) * inv(fact[n-k]) 로 구할 수 있다. 전처리는 O(n), 쿼리 한번은 모듈러 역원을 두번 계산하면 되니까 O(logP). 즉 q개의 쿼리를 O(n + qlogP)에 계산 가능.
전처리 단계에서 fact_inv[i] = inv(i!) 도 모두 계산해 둔다면, 쿼리 한번을 O(1)에 계산할수 있다.
fact_inv[i] 테이블을 채우는 첫번째 방법은, inv( (n-1)!) = n*inv(n!) 임을 이용하는 방법이다. inv(n!)를 O(logP)에 구하면 나머지 셀들 (fact_inv[1] ~ fact_inv[n-1])을 채우는 것은 O(n)에 가능하다. 그러면 총 전처리 시간은 O(logP + n)
다른 방법은 inv(n) = -(p / i) * inv(p % i)을 이용하는 방법이다. inv[i] 테이블을 O(n)에 채우고, fact_inv[i] = fact_inv[i-1] * inv[i] 를 이용해서, fact_inv테이블 역시 O(n)에 채울 수 있다. 앞서 말했듯 logP가 n에 비해서 작으니, 윗 방법이나 이 방법이나 시간복잡도 차이는 없다고 봐도 되고, 코딩이 편한 쪽을 쓰면 된다
이렇게 fact_inv 테이블 역시 O(n)에 채워 둔다면, q개의 쿼리를 O(n+q)에 처리 가능하다.
코드
관련 문제
특수한 경우들 - 모듈러스 M이 n보다 큰 소수가 아닌 경우
아예 이것을 물어보는 정수론 문제가 아니고서, 일반적인 문제에서는 잘 쓰일 일이 없는 경우들이다.
정확히는 모듈러스 M이 k!과 서로소가 아닐 때에는 inv(k!)를 구하는 것이 불가능해서 위의 방법을 사용할 수 없어서, 아래에 설명된 방법들을 써야 한다.
분자와 분모를 모두 소인수분해해서 계산하는 방법
코드
for p in numtheory.prime_list(N):
prime_pow = p
e = 0
while prime_pow <= N:
e += (N // prime_pow) - (K // prime_pow) - ((N - K) // prime_pow)
prime_pow *= p
ans = ans * pow(p, e, M) % M
관련 문제
모듈러스가 n보다 작은 소수일 때
정확히는 k보다 작은 소수일 때. M이 n보다 작더라도 k보다 크면 k!의 역원이 존재하므로 그냥 계산 가능
두가지 방법이 있다.
뤼카의 정리를 이용하는 방법과 n!_P % P 를 구해서 사용하는 것. 뤼카의 정리를 이용하는 방법이 시간도 효율적이고 코드도 간단하다.
-
C(n,k) % P 를, $log_P{n}$ 개의 C(ni, ki)의 곱으로 계산할 수 있다. (ni < P).
이 방식을 이용하면, O(P + logn/logP)로 계산 가능. 쿼리가 q개가 들어온다면, O(P + q*logn/logP)에 처리 가능.
2. n!_P % m 을 구해서 이용하는 방법. (n!_P = (n!)/(P^c(n)). n!에서 P성분을 모두 제거한 것)
n!/(k!(n-k)!) = n!_P / (k!_P)( (n-k)!_P) * P^(c(n)-c(k)-c(n-k)) 로 고쳐써서 계산하는 방법이다. 이렇게 바꿔쓰면, k!_P와 (n-k)!_P는 P와 서로소이므로 P에 대한 모듈러 역원이 존재하므로 계산이 가능하다.
n!_P 는
팩토리얼에서 설명했던 방법으로 O(P)의 전처리이후, O(logn/logP)에 구할 수 있다. 이렇게 n!_P, k!_p, (n-k)!_P 를 모두 O(P + logn/logP)에 구한다.
역시
팩토리얼애서 설명한 방법으로 c(n), c(k), c(n-k)를 구한다. 시간복잡도는 O(logn/logP),
이제 c(n)-c(k)-c(n-k) > 0 이면 n!/(k!(n-k)!) 이 P로 나눠져 떨어진다는 의미이므로 답은 그냥 0. c(n)-c(k)-c(n-k) = 0일때만 계산하면 되고, 이 경우는 n!/(k!(n-k)!) = n!_P / (k!_P)( (n-k)!_P) 이다.
k!_P와 (n-k)!_P의 모듈러 역원을 구하는 데에 O(logP), 전체 시간복잡도는 O(P + logn/logP + logP)
쿼리가 여러개 들어온다면 O(P + q*(logn/logP + logP)).
1번이 좀더 코드도 짧고, 속도도 1번이 좀 더 빠르다.
코드
from teflib import combinatorics
comb_table = combinatorics.CombTable(M - 1, M)
answer = 1
while N > 0:
N, n_mod = divmod(N, M)
K, k_mod = divmod(K, M)
answer *= comb_table.get(n_mod, k_mod)
answer %= M
관련 문제
모듈러스 M이 p^e 형태일때
코드
fact, inv = precompute_for_prime_pow(p, e):
def comb_mod_prime_pow(n, k, p, e):
"""Returns C(n, k) % (p^e)."""
mod = p**e
r = n - k
answer = 1
t = 0
while n > 0:
answer *= fact[n % mod] * inv[k % mod] * inv[r % mod]
if (n // mod + k // mod + r // mod) % 2 == 1:
answer *= -1
n //= p
k //= p
r //= p
t += n - k - r
return answer * pow(p, t) % mod
모듈러스가 임의의 합성수일 때
이제 드디어 마지막 단계.
M을 p1^e1 * p2^e2 * … 형태로 소인수분해하고, C(n,k) % p1^e1, C(n,k) % p2^e2, …을 각각 계산하고
중국인의 나머지 정리를 적용해서 최종 답을 구한다.
C(n,k) % p^e 을 구하는 방법은 위에서 설명한대로 뤼카의 정리나 (N!_p % p^e)을 이용해서 구하면 된다. e=1일때는 뤼카의 정리를 쓰고, e>1일때는 (N!_p % p^e) 방법으로 구하는 것이 좀 더 빠르긴 한데, 두가지 방법을 다 구현하는게 귀찮으면 그냥 e에 관계 없이 (N!_p % p^e) 방법으로 풀어도 된다.
코드
def precompute(MOD_FACTORS):
facts = {}
invs = {}
for p, e in MOD_FACTORS.items():
mod = p**e
facts[mod], invs[mod] = precompute_for_prime_pow(p, e)
return facts, invs
def comb_mod_m(N, K, MOD_FACTORS):
rems = []
mods = []
for p, e in MOD_FACTORS.items():
mod = p**e
rems.append(comb_mod_prime_pow(N, K, p, e, facts[mod], invs[mod]))
mods.append(mod)
print(numtheory.linear_congruences(rems, mods))
관련 문제