2

当n，r，m是非常大的数字时，如何在编程中计算nCr % m（即“（n选择r）模数m”）？

4

4 回答 4

1

我假设您想计算 nCr 类型的大量数据并搜索一些可以优化的库函数。在这种情况下，您可以查看GNU Scientific Library。

于 2012-11-08T06:53:28.127 回答
1

首先，这取决于您所说的“非常大”是什么意思。如果您的问题仅仅是对于标准 64 位整数来说中间值变得太大，那么您可以使用gmp或BigInteger之类的东西。

所涉及的数字可能会变得如此之大，以至于您将耗尽内存或耐心，并且您将无法计算中间值以完成精度。在这种情况下，最好的办法是首先确定每个阶乘的素因式分解，使用这些因式分解来确定二项式的因式分解，然后在每个中间步骤使用模数乘以这个因式分解。

您将需要一个不超过 n 的所有素数的列表。

以下是伪代码。我正在使用int，但您应该将其替换为您正在使用的任何大型库。

int factorial_prime_power(int f, int p) {
    // When p is prime, returns k such that p^k divides f!
    int k = 0;
    while (f > p) {
        f = f / p;
        k = k + f;
    }
    return k;
}
int binomial_prime_power(int n, int r, int p) {
    // when p is prime, returns k such that p^k divides nCr
    return factorial_prime_power(n,p) - factorial_prime_power(r,p) - factorial_prime_power(n-r,p);
}
int powmod(int p, int k, int m) {
    // quickly calculates p^k mod m
    int res = 1;
    int q = p;
    int j = k;
    while (j > 0) {
        // invariant:  p^k is congruent to res * q^j
        if (j is odd) {
            res = (res * q) % m;
            j = (j-1)/2;
        } else {
            j = j / 2;
        }
        q = (q * q) % m;
    }
    return res;
}
int big_binomial(int n, int r, int m) {
    if (n < r or r < 0) {
        return 0;
    }
    int res = 1;
    for(p in all primes from 2 to n) {
        k = binomial_prime_power(n,r,p);
        res = (res * powmod(p,k,m)) % m;
    }
    return res;
}
于 2012-11-11T15:25:38.380 回答
1

有一些技术可以降低计算以整数为模的二项式系数的复杂性m，前提是m的主要因素都不是太大。m以过高的功率划分。

第一步是因式分解m，

m = ∏ (p_i ^ e_i)

然后计算这些素数中的每一个的模的二项式系数，并将结果与​​中国剩余定理相结合。

以素数 ( ) 为模的二项式系数的计算e_i == 1可以比一般情况更容易计算，参见。例如这个答案，但它也可以包含在那里。

对于素数p和n >= 0，让我们定义

F(p, n) =   ∏ k   = p^(n/p) * (n/p)!
         1<=k<=n
           p | k

和

G(p, n) =    ∏ k
          1<=k<=n
         gcd(k,n)=1

然后我们有

n! = F(p, n) * G(p, n)

并且迭代地，对(n/p)!出现在 中使用相同的拆分F(p, n)，

           m
n! = p^K * ∏ G(p, n/(p^j))
          j=0

哪里p^m <= n < p^(m+1)。中的所有因子G(p, x)都与 互质p^e，因此二项式系数的分母中的相应因子可以取模p^e，如果我们找到一种有效的G(p, x)模数计算方法p^e，我们就有一种有效的方法来计算二项式系数的模数p^e。

对于二项式系数，我们有

n! / (r! * (n-r)!) = p^M * (∏ G(p, (n/p^j)) * [ ∏ G(p, r/(p^j)) * ∏ G(p, (n-r)/(p^j)) ]^(-1)

让H(p, e, n) = G(p, n) % (p^e). 关键是所有数的互质p^e不超过的乘积p^e相当简单。它与-1模一致p^e，除非p = 2和e > 2，在这种情况下它与 1 一致。

所以

H(p, e, n) ≡ (-1)^(n/(p^e)) * H(p, e, n % (p^e)) (mod p^e)

（除非p = 2和e > 2，在这种情况下第一个因子是 1），我们只需要计算H(p, e, k)，0 <= k < p^e然后我们可以查找结果。

代码：

// invert k modulo p, k and p are supposed coprime
unsigned long long invertMod(unsigned long long k, unsigned long long p) {
    unsigned long long q, pn = 1, po = 0, r = p, s = k;
    unsigned odd = 1;
    do {
        q = r/s;
        q = pn*q + po;
        po = pn;
        pn = q;
        q = r%s;
        r = s;
        s = q;
        odd ^= 1;
    }while(pn < p);
    return odd ? p-po : po;
}

// Calculate the binomial coefficient (n choose k) modulo (prime^exponent)
// requires prime to be a prime, exponent > 0, and 0 <= k <= n,
// furthermore supposes prime^exponent < 2^32, otherwise intermediate
// computations could have mathematical results out of range.
// If k or (n-k) is small, a direct computation would be more efficient.
unsigned long long binmod(unsigned long long prime, unsigned exponent,
                          unsigned long long n, unsigned long long k) {
    // The modulus, prime^exponent
    unsigned long long ppow = 1;
    // We suppose exponent is small, so that exponentiation by repeated
    // squaring wouldn't gain much.
    for(unsigned i = 0; i < exponent; ++i) {
        ppow *= prime;
    }
    // array of remainders of products
    unsigned long long *remainders = malloc(ppow * sizeof *remainders);
    if (!remainders) {
        fprintf(stderr, "Allocation failure\n");
        exit(EXIT_FAILURE);
    }
    for(unsigned long long i = 1; i < ppow; ++i) {
        remainders[i] = i;
    }
    for(unsigned long long i = 0; i < ppow; i += prime) {
        remainders[i] = 1;
    }
    for(unsigned long long i = 2; i < ppow; ++i) {
        remainders[i] *= remainders[i-1];
        remainders[i] %= ppow;
    }

    // Now to business.
    unsigned long long pmult = 0, ntemp = n, ktemp = k, mtemp = n-k,
                       numer = 1, denom = 1, q, r, f;
    if (prime == 2 && exponent > 2) {
        f = 0;
    } else {
        f = 1;
    }
    while(ntemp) {
        r = ntemp % ppow;
        q = ntemp / ppow;
        numer *= remainders[r];
        numer %= ppow;
        if (q & f) {
            numer = ppow - numer;
        }
        ntemp /= prime;
        pmult += ntemp;
    }
    while(ktemp) {
        r = ktemp % ppow;
        q = ktemp / ppow;
        denom *= remainders[r];
        denom %= ppow;
        if (q & f) {
            denom = ppow - denom;
        }
        ktemp /= prime;
        pmult -= ktemp;
    }
    while(mtemp) {
        r = mtemp % ppow;
        q = mtemp / ppow;
        denom *= remainders[r];
        denom %= ppow;
        if (q & f) {
            denom = ppow - denom;
        }
        mtemp /= prime;
        pmult -= mtemp;
    }
    // free memory before returning, we don't use it anymore
    free(remainders);
    if (pmult >= exponent) {
        return 0;
    }
    while(pmult > 0) {
        numer = (numer * prime) % ppow;
        --pmult;
    }
    return (numer * invertMod(denom, ppow)) % ppow;
}

它n choose k modulo p^e分O(p^e + log n)步计算。

于 2012-11-08T20:10:47.957 回答
0

分块做。

例如：n!/(r!(nr)!) = 1*2*...*n/(1*2*...r*1*2*...*(nr))=1/ 1*1 * 2/2*2 * 3/3*3 * ... *... 。

每个块都很容易计算，因此您应该避免计算中的溢出。计算块并将您当前的结果乘以它。

此外，在此之前对 n 和 r 进行取模是值得的。

于 2012-11-08T07:01:02.147 回答