当n非常大时，高效计算nCr mod p。

Question

我需要高效计算nCr mod p。目前，我已经编写了这段代码，但它超出了时间限制。请建议更优的解决方案。

对于我的情况，p = 10^9 + 7，1 ≤ n ≤ 100000000 我还要确保没有溢出，因为nCr mod p保证适合32位整数，但是n!可能会超过限制。

def nCr(n,k):
    r = min(n-k,k)
    k = max(n-k,k)
    res = 1
    mod = 10**9 + 7

    for i in range(k+1,n+1):
        res = res * i
        if res > mod:
            res = res % mod

    res = res % mod
    for i in range(1,r+1):
        res = res/i
    return res

PS：我认为我的代码可能不完全正确。然而，它似乎能够正确地处理小n。如果有错误，请指出!

- OneMoreError

---

你使用的是哪个版本的Python？ - inspectorG4dget

我正在使用Python 2.7.2。 - OneMoreError

你为什么担心溢出？Python 整数类型没有固定的存储空间；它将使用所需的存储空间。 - Ramchandra Apte

@RamchandraApte：虽然那是真的，但在Python中仍然有可能导致整数溢出。 - inspectorG4dget

@inspectorG4dget 嗯，我只知道在Python 3中使用某些C函数，例如math模块中的函数。对于Python 2不太确定。 - Ramchandra Apte

Answer 1

来自http://apps.topcoder.com/wiki/display/tc/SRM+467：

long modPow(long a, long x, long p) {
    //calculates a^x mod p in logarithmic time.
    long res = 1;
    while(x > 0) {
        if( x % 2 != 0) {
            res = (res * a) % p;
        }
        a = (a * a) % p;
        x /= 2;
    }
    return res;
}

long modInverse(long a, long p) {
    //calculates the modular multiplicative of a mod m.
    //(assuming p is prime).
    return modPow(a, p-2, p);
}
long modBinomial(long n, long k, long p) {
// calculates C(n,k) mod p (assuming p is prime).

    long numerator = 1; // n * (n-1) * ... * (n-k+1)
    for (int i=0; i<k; i++) {
        numerator = (numerator * (n-i) ) % p;
    }

    long denominator = 1; // k!
    for (int i=1; i<=k; i++) {
        denominator = (denominator * i) % p;
    }

    // numerator / denominator mod p.
    return ( numerator* modInverse(denominator,p) ) % p;
}

注意我们使用 modpow(a, p-2, p) 来计算模逆。这是根据费马小定理，即 (a^(p-1) ≡ 1 (mod p))，其中p是质数。因此，它意味着 (a^(p-2) ≡ a^(-1) (mod p))。

C++ 转 Python 应该很容易 :)

- batbaatar

---

另外，请注意在Python中，pow()已经以modPow的形式提供。 - Ramchandra Apte

帮了我很大的忙。在其他地方很难找到这个实现方式。 - ryan1234

Answer 2

关于最后一个问题：我认为你代码中的错误在于计算乘积，对k取模，然后再将结果除以r!。这与在对k取模之前进行除法不同。例如，3*4 / 2 (mod 10) != 3*4 (mod 10) / 2。

- Armin Rigo

Answer 3

从 sympy 1.13 开始（截至今天尚未发布），可以使用函数 binomial_mod 来实现这一功能。这是合并的 PR：https://github.com/sympy/sympy/pull/24891 它使用了广义的 Lucas 定理。

- Rexcirus