技術メモ

役に立てる技術的な何か、時々自分用の覚書。幅広く色々なことに興味があります。

nCk (mod m) の求め方 [n,kが凄く大きい場合]

Python 数学

大きいｎに対してコンビネーションを求める

凄く大きいnに対してCombinationを求める時、例えば {}_{1000} \mathrm{C}_{400}を求めようと思っても愚直に計算するにはプログラムで扱える桁数を超えてしまう。
そこで、mで割った余り(modular) {}_n \mathrm{C}_k \pmod{m}を求めさせることがある。

pythonで書いたのでサンプルコードを載せておく。
ただし、これを使うにはmはnを超える素数でないといけない。

サンプルコード

再帰を使っている。

def mod_combination(n, k, mod):
    # nCk (mod m)
    def mod_permutation(n, k, mod):
        if n<=k:
            return 1
        else:
            return (n * mod_permutation(n-1,k,mod))%mod

    def mod_inv_permutation(k, mod):
        k, mod = int(k), int(mod)
        if k<=1:
            return 1
        else:
            return (pow(k,mod-2,mod) * mod_inv_permutation(k-1, mod))%mod

    return (mod_permutation(n,n-k,mod) * mod_inv_permutation(k, mod))%mod

理屈

割り算に対してモジュラ（余り）を計算するにはどうするかという問題がある。
それにはモジュラ逆数を考える必要がある。
プログラムでは

pow(k, mod-2, mod)

にあたる部分でモジュラー逆数を使っている。

モジュラ逆数とはa^{-1}\equiv x \pmod{m}となるxのこと。

結論から言うと

素数 mと互いに素な aに対して
a^{-1}\equiv a^{m-2} \pmod{m}

となることを使っている。

(証明)
オイラー関数 \phi(m)を用いて、mが素数の時
\displaystyle \left\{ \begin{array}{1} a^{\phi(m)}\equiv 1 \pmod{m}\\ \phi(m) = m-1 \end{array} \right .
なので
a^{\phi(m) - 1} \equiv a^{-1} \pmod{m} \\ a^{m-2} \equiv a^{-1} \pmod{m}
(Q.E.D)

したがって {}_n \mathrm{C}_k (n < m)に対して
\begin{eqnarray} {}_n \mathrm{C}_k & = & n(n-1)\dots(n-k) \times k^{-1}(k-1)^{-1}\dots 1^{-1} \\ & \equiv & n(n-1)\dots(n-k)\times k^{m-2}(k-1)^{m-2}\dots 1^{m-2} \pmod{m} \end{eqnarray}


上のプログラムでは

def mod_permutation(n, k, mod)

で n(n-1)\dots(n-k) \pmod{m} 　を

def mod_inv_permutation(n, k, mod)

で k^{m-2}(k-1)^{m-2}\dots 1^{m-2} \pmod{m}　を
計算している。

おまけ

素数を確認したい時pythonではsympyを使おう

from sympy import isprime
print(isprime(1000000007)) #True

参考サイト

How to find the value of nCr%1000000007 for the large number n (n <= 10^6) in C - Quora
モジュラ逆数 - Wikipedia

mが素数以外の場合も :
nCr mod mの求め方 - uwicoder - アットウィキ