• 정의와 설명은 ko:이항 계수 참고.
  • $ {\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일 때는 설명이 다르다
    • 한국어 위키에서는 k<0 일때에도 C(n,k)=0이라고 되어있다. 영문 위키나 MathWorld에서는 n과 k가 음이 아닌 정수일때만 C(n,k)가 정의된다고 한다.
    • python의 math.comb 함수는 n이나 k가 음수이면 ValueError를 발생시키고, k>n일때는 0을 리턴한다.
• 여기저기에서 너무도 많이 사용된다. 이항 계수를 포함하는 형태로 정의되는 다른 수열들도 많다. ex) 카탈랑 수 (Catalan number), 스털링 수
• C(n,k)의 값은 빠르게 커지므로, n>=35면 32비트 정수형의 범위를 넘어가고 n>=68이면 64비트 정수형의 범위를 넘어가게 된다. 그래서 보통 문제에서는 모듈러를 취한 값을 요구하게 된다.
• C(n,k)에서 k>n-k 인 경우에는 C(n,n-k)을 구해도 값이 똑같으므로, 아래부터는 항상 k<n-k임을 가정하고 설명하겠다.

• 쿠머의 정리
  • Kummer's theorem, https://jjycjnmath.tistory.com/532
  • 어떤 소수 p에 대해서, p^k가 C(n,m)을 나누어 떨어지게 하는 가장 큰 k값을 찾는 문제.
    • 표현이 긴데, 원래는 대응되는 용어가 있다. 위에서 말한 k를 C(n,m)의 p-adic valuation 이라고 부르고 $\nu_{p}({\binom {n}{m}}) $ 이라고 쓴다
      • $ \nu_{p}(n) = \begin{cases} \max \{\nu \in N \, : \, p^{\nu} \mid n \} & \text{if} \;\; n \neq 0 \\ \infty & \text{if} \;\; n = 0 \end{cases} $
  • $ \nu _{p}\left({\binom {n}{m}}\right)={\dfrac {S_{p}(m)+S_{p}(n-m)-S_{p}(n)}{p-1}} $
    • 간단히 말하면 $\nu_{p}({\binom {n}{m}}) $ 은 p진법에서 n-m에 m 을 더할 때 생기는 받아올림의 갯수와 같다.
  • 직관적인 활용예로, C(n,m)의 홀짝성을 찾는 문제가 있다.
    • p=2일때 $\nu_{p}({\binom {n}{m}}) $ 이 0이면 홀수, 1이상이면 짝수라는 뜻이므로, 2진수끼리 더할때 받아올림이 생기는지 없는지만 확인하면 된다. 받아올림이 생기는 경우는 같은 자리에 bit값이 둘다 1인 경우뿐이므로, 결국 이것은 그냥 (n-m)과 m의 bitwise and 연산 결과가 0인지 아닌지만 보면 된다. 0이면 홀수, 1이면 짝수이다.
• 뤼카의 정리
  • Lucas's theorem
  • 사실상 쿠머의 정리를 조금 다른 형태로 나타낸 정리이다.
  • $ {\binom {m}{n}}\equiv \prod _{i=0}^{k}{\binom {m_{i}}{n_{i}}}{\pmod {p}} $
    • m_i와 n_i 는 m과 n을 p진법으로 나타냈을때 i번째 자리의 숫자이다.
  • 이항계수를 소수로 나눈 나머지를 구해야 할때, 더 작은 값에 대한 이항계수들의 곱으로 구할수 있는 방법을 제공한다. 이는 실제로 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를 사용해야 한다.

• 이항계수를 삼각형 모양으로 펼쳐놓은 것이 ko:파스칼의 삼각형이고, 여기서의 핵심은 C(n,k) = C(n-1,k-1) + C(n-1,k) 의 점화식이다.
• 이 방식의 특징은 덧셈만을 갖고서 이항계수를 구할 수 있다는 것.
  • 큰값을 구하기 위해서 BigInteger를 직접 구현해서 사용해야 할때, 이 방법으로는 BigInteger의 덧셈 연산만 구현해도 적용 가능하다.
  • 모듈러를 취한 값을 구해야 할때도, 모듈러 나눗셈을 사용하지 않고도 계산이 가능하다. 따라서 모듈러스의 소수 여부나 모듈러스의 크기에 관계 없이 사용이 가능하다
• 하지만 DP를 사용해도 시간복잡도가 O(n^2)으로 비효율적이라서 실제로 쓸 일은 거의 없다.
  • 특히 파이썬에서는 이 방법으로 돌릴 수 있을 정도의 n의 범위에서는, 모듈러를 취한 값을 계산해야 할 경우조차도, 그냥 정확한 값을 math.comb를 써서 구하고 마지막에 모듈러를 취해버리는 것이 더 빠르다.

• 일반적으로 많이 쓰게 될 유형들이다. P>n 인 소수 P에 대해서 C(n,k) % P 를 구하는 방법이다.
• 아래 사이트들에 있는 내용들을 재정리한 것이다.
  • https://rkm0959.tistory.com/182
  • https://koosaga.com/63
  • https://cp-algorithms.com/combinatorics/binomial-coefficients.html

• 모듈러 나눗셈을 이용해서, 팩토리얼 형태의 공식대로 계산하는게 가장 깔끔하다. 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 정도까지도 위 방법에 비해서 속도가 크게 차이나지 않는다.
  • python 3.11에는 math.comb의 속도를 더욱 빠르게 하는 업데이트가 적용될 예정이다 (https://bugs.python.org/issue37295)

teflib.combinatorics.comb

• 이항 계수 3

• 그냥 C(n,k)를 한번만 구하고 끝인 것이 아니라, C(n1,k1), C(n2,k2), …, C(nq,kq) 를 모두 구해야 하는 경우이다 (ni ≤ n)
• 단순히 값을 한번만 계산할 때의 방법을 q번 반복하면 O(q(n+logP)). 이건 너무 비효율적d이다.
• 전처리 단계에서 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)에 처리 가능하다.

teflib.combinatorics.CombTable

• 이항 계수와 쿼리

• 아예 이것을 물어보는 정수론 문제가 아니고서, 일반적인 문제에서는 잘 쓰일 일이 없는 경우들이다.
• 정확히는 모듈러스 M이 k!과 서로소가 아닐 때에는 inv(k!)를 구하는 것이 불가능해서 위의 방법을 사용할 수 없어서, 아래에 설명된 방법들을 써야 한다.

• 모듈러스에 관계 없이 사용할 수 있는 방법으로 n이 적당히 작을때 유리하다. 이후에 설명할 방법들과는 달리 시간복잡도가 M이 아닌 n에 관한 식으로 나온다.
• 소인수분해 방법을 이용하여, n!, k!, (n-k)!을 모두 소인수분해 한다. 그리고 소인수분해 결과를 이용해서 n!/k!(n-k)! 을 소인수분해 형태, p1^e1*p2^e2*…*pn^en 의 형태로 계산한다. 그 다음 이것을 M으로 나눈 나머지를 (p1^e1 % M) * (p2^e2 % M) * (p3^e3 % M) * … 으로 계산한다.
  • O(nloglogn)의 전처리를 통해 n 이하의 모든 소수를 구해둔다면, 각각을 소인수 분해하는 것은 O(n)이하에 가능하다. 소인수 분해 결과를 이용해서 n!/k!(n-k)!를 구하는 것, 그것을 M으로 나눈 나머지를 진짜로 계산하는 것은 모두 O(n)이하이므로, 시간 복잡도는 O(nloglogn)에 바운드된다

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

• 이항 계수 5

• 정확히는 k보다 작은 소수일 때. M이 n보다 작더라도 k보다 크면 k!의 역원이 존재하므로 그냥 계산 가능
• 두가지 방법이 있다. 뤼카의 정리를 이용하는 방법과 n!_P % P 를 구해서 사용하는 것. 뤼카의 정리를 이용하는 방법이 시간도 효율적이고 코드도 간단하다.
• 1. 뤼카의 정리를 이용하는 방법.
  • C(n,k) % P 를, $log_P{n}$ 개의 C(ni, ki)의 곱으로 계산할 수 있다. (ni < P).
    • 공식을 적용한 코드는 짜기 어렵지 않다. 여러번의 C(ni, ki) % P를 빠르게 계산하기 위해서, 값을 여러번 계산해야 할 때에서 설명한대로, P크기의 fact과 inv 테이블을 미리 만들어 놓고 사용하는 게 좋다.
    • 다만 공식대로 ni과 ki를 구하면, ni<ki 인 경우도 존재할 수 있다. 이 경우에 C(ni, ki)가 0이 되도록 처리하는 것에 주의할 것.
  • 이 방식을 이용하면, 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)).
    • 만약 n!_P와 모듈러 역원들, 그리고 c(n)을 모두 다 O(n)에 전처리해둔다면 각 쿼리는 O(1)에 계산 가능하므로 O(n + q)로 바뀐다. 그러나 logn/logP + logP 이 거의 상수에 가깝기 때문에 그냥 계산해도 거의 O(P+q)이므로 굳이 이렇게 할 필요는 없을듯.
• 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

• 이항 계수 4

• 여기서 기본 조건은 당연히 p<k일때이다. p>n이면 p^e과 n!이 서로소이므로 모듈러스가 n보다 큰 소수일때와 마찬가지로 그냥 모듈러 역원을 구해서 계산하면 된다.
• 모듈러스가 n보다 작은 소수일 때에서의 두가지 방법을 각각 확장해서 적용 가능하다.
• 1. 일반화된 뤼카의 법칙을 사용하는 방법이다. 다만 이 방법은 위키피디아나 cp-algorithms에서도 설명 없이 논문 링크만 걸려있어서 접근성이 좀 떨어졌는데, 최근에 방법을 설명한 글이 https://velog.io/@stripe2933/bj14854a에 올라왔다. (그러나 자세히 읽어보지는 않았다..ㅜ)
• 2. 위에서와 마찬가지로 n!/(k!(n-k)!) = n!_P / (k!_P)( (n-k)!_P) * P^(c(n)-c(k)-c(n-k)) 으로 변환해서 푸는 방법이다.
  • 다만 위에서는 n!_P % P를 계산해야 했고, 여기에서는 n!_P % (p^e)를 계산해야 하는 차이가 있다. 팩토리얼의 방법으로 O(p^e)=O(M)의 전처리 이후, O(logn/logP)의 계산이 필요하다.
  • c(n),c(k),c(n-k)를 구하는 것은 똑같다. O(logn/logP). 그리고 역시 c(n)-c(k)-c(n-k) ≥ e 이면, 답은 0이다.
  • k!_P와 (n-k)!_P의 모듈러 역원을 구하는 데에 O(logM), P^(c(n)-c(k)-c(n-k)를 구하는 것은 O(loge), 전체 시간복잡도는 O(M + logn/logP + logM + loge) 정도로 나온다. 사실 M을 제외한 나머지 텀들은 거의 상수에 가깝다..
  • 쿼리가 여러개 들어온다면 O(M + q*(logn/logP + logM + loge)). 만약 n!_P와 모듈러 역원들, 그리고 c(n)을 모두 다 O(n)에 전처리해둔다면 각 쿼리는 O(1)에 계산 가능하므로 O(n + q)로 바꿀수 있지만, 그냥 계산하는 것도 거의 O(M+q)에 가까우므로 비효율적.

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))

• 돌아온 떡파이어
  • M=p1*p2의 형태. 뤼카의 정리와 중국인의 나머지 정리만으로 풀수 있다
• 이항 계수 6
  • M=p1^e1*p2*p3의 형태. 이제껏 나온 모든 방법을 사용해서 푼다.