47

我所说的“大 n”是指数以百万计的东西。p 是素数。

我已经尝试过 http://apps.topcoder.com/wiki/display/tc/SRM+467 但该功能似乎不正确(我用 144 选择 6 mod 5 对其进行了测试,当它应该给我时它给了我 0 2)

我试过 http://online-judge.uva.es/board/viewtopic.php?f=22&t=42690 但我不完全理解

我还制作了一个使用逻辑 (combinations(n-1, k-1, p)%p + combination(n-1, k, p)%p) 的记忆递归函数,但它给了我堆栈溢出问题,因为n 很大

我试过卢卡斯定理,但它似乎很慢或不准确。

我要做的就是创建一个快速/准确的 n 为大 n 选择 k mod p。如果有人可以帮助我展示一个很好的实现,我将不胜感激。谢谢。

根据要求,对于大 n 命中堆栈溢出的记忆版本:

std::map<std::pair<long long, long long>, long long> memo;

long long combinations(long long n, long long k, long long p){
   if (n  < k) return 0;
   if (0 == n) return 0;
   if (0 == k) return 1;
   if (n == k) return 1;
   if (1 == k) return n;

   map<std::pair<long long, long long>, long long>::iterator it;

   if((it = memo.find(std::make_pair(n, k))) != memo.end()) {
        return it->second;
   }
   else
   {
        long long value = (combinations(n-1, k-1,p)%p + combinations(n-1, k,p)%p)%p;
        memo.insert(std::make_pair(std::make_pair(n, k), value));
        return value;
   }  
}
4

3 回答 3

59

因此,这是解决问题的方法。

你当然知道公式:

comb(n,k) = n!/(k!*(n-k)!) = (n*(n-1)*...(n-k+1))/k! 

(参见http://en.wikipedia.org/wiki/Binomial_coefficient#Computing_the_value_of_binomial_coefficients

你知道如何计算分子:

long long res = 1;
for (long long i = n; i > n- k; --i) {
  res = (res * i) % p;
}

现在,由于 p 是素数,因此与 p 互质的每个整数的倒数都已明确定义,即可以找到-1 。这可以使用费马定理 a p-1 =1(mod p) => a*a p-2 =1(mod p) 等 a -1 =a p-2 来完成。现在您需要做的就是实现快速求幂(例如使用二进制方法):

long long degree(long long a, long long k, long long p) {
  long long res = 1;
  long long cur = a;

  while (k) {
    if (k % 2) {
      res = (res * cur) % p;
    }
    k /= 2;
    cur = (cur * cur) % p;
  }
  return res;
}

现在您可以将分母添加到我们的结果中:

long long res = 1;
for (long long i = 1; i <= k; ++i) {
  res = (res * degree(i, p- 2)) % p;
}

请注意,我在任何地方都使用 long long 以避免类型溢出。当然你不需要做k幂运算——你可以计算 k!(mod p) 然后只除一次:

long long denom = 1;
for (long long i = 1; i <= k; ++i) {
  denom = (denom * i) % p;
}
res = (res * degree(denom, p- 2)) % p;

编辑:根据@dbaupp 的评论 if k >= p the k! 将等于 0 模 p 并且 (k!)^-1 不会被定义。为了避免这种情况,首先计算 p 在 n*(n-1)...(n-k+1) 和 k 中的程度!并比较它们:

int get_degree(long long n, long long p) { // returns the degree with which p is in n!
  int degree_num = 0;
  long long u = p;
  long long temp = n;

  while (u <= temp) {
    degree_num += temp / u;
    u *= p;
  }
  return degree_num;
}

long long combinations(int n, int k, long long p) {
  int num_degree = get_degree(n, p) - get_degree(n - k, p);
  int den_degree = get_degree(k, p);

  if (num_degree > den_degree) {
    return 0;
  }
  long long res = 1;
  for (long long i = n; i > n - k; --i) {
    long long ti = i;
    while(ti % p == 0) {
      ti /= p;
    }
    res = (res * ti) % p;
  }
  for (long long i = 1; i <= k; ++i) {
    long long ti = i;
    while(ti % p == 0) {
      ti /= p;
    }
    res = (res * degree(ti, p-2, p)) % p;
  }
  return res;
}

编辑:还有一个可以添加到上述解决方案的优化 - 我们可以计算 k!(mod p) 然后计算该数字的倒数,而不是计算 k! 中每个倍数的倒数。因此,我们只需要支付一次取幂的对数。当然,我们必须再次丢弃每个倍数的 p 除数。我们只需要改变最后一个循环:

long long denom = 1;
for (long long i = 1; i <= k; ++i) {
  long long ti = i;
  while(ti % p == 0) {
    ti /= p;
  }
  denom = (denom * ti) % p;
}
res = (res * degree(denom, p-2, p)) % p;
于 2012-04-12T06:17:03.733 回答
13

对于大k,我们可以通过利用两个基本事实来显着减少工作:

  1. 如果p是素数,则p的素数因式分解中的指数n!由 给出(n - s_p(n)) / (p-1),其中是基本表示s_p(n)中的数字之和(因此对于,它是 popcount)。因此,在 的素数分解中的指数是,特别是,当且仅当在基数中执行加法时没有进位(指数是进位的数量)时,它才为零。npp = 2pchoose(n,k)(s_p(k) + s_p(n-k) - s_p(n)) / (p-1)k + (n-k)p

  2. 威尔逊定理: p是素数,当且仅当(p-1)! ≡ (-1) (mod p)

p因式分解中的指数n!通常由下式计算

long long factorial_exponent(long long n, long long p)
{
    long long ex = 0;
    do
    {
        n /= p;
        ex += n;
    }while(n > 0);
    return ex;
}

choose(n,k)对by的可分性检查p并不是绝对必要的,但首先进行检查是合理的,因为通常是这样,然后工作量就会减少:

long long choose_mod(long long n, long long k, long long p)
{
    // We deal with the trivial cases first
    if (k < 0 || n < k) return 0;
    if (k == 0 || k == n) return 1;
    // Now check whether choose(n,k) is divisible by p
    if (factorial_exponent(n) > factorial_exponent(k) + factorial_exponent(n-k)) return 0;
    // If it's not divisible, do the generic work
    return choose_mod_one(n,k,p);
}

现在让我们仔细看看n!。我们将数字≤ n分成 的倍数p和互质的数字p。和

n = q*p + r, 0 ≤ r < p

p贡献的倍数p^q * q!。这些数字互质以p贡献(j*p + k), 1 ≤ k < pfor0 ≤ j < q的乘积和 的乘积(q*p + k), 1 ≤ k ≤ r

对于互质的数字,p我们只会对贡献模感兴趣p。每个完整的运行j*p + k, 1 ≤ k < p都与(p-1)!模数一致p,因此它们总共产生(-1)^q模数的贡献p。最后一次(可能)不完整的运行产生r!模数p

所以如果我们写

n   = a*p + A
k   = b*p + B
n-k = c*p + C

我们得到

choose(n,k) = p^a * a!/ (p^b * b! * p^c * c!) * cop(a,A) / (cop(b,B) * cop(c,C))

其中是cop(m,r)与 互质的所有数的乘积。p≤ m*p + r

有两种可能性,a = b + cA = B + C,或a = b + c + 1A = B + C - p

在我们的计算中,我们已经预先排除了第二种可能性,但这不是必须的。

在第一种情况下,p取消的明确权力,我们留下

choose(n,k) = a! / (b! * c!) * cop(a,A) / (cop(b,B) * cop(c,C))
            = choose(a,b) * cop(a,A) / (cop(b,B) * cop(c,C))

任何p除法choose(n,k)都来自choose(a,b)- 在我们的例子中,不会有,因为我们之前已经消除了这些情况 - 并且,虽然cop(a,A) / (cop(b,B) * cop(c,C))不必是整数(考虑 eg choose(19,9) (mod 5)),当考虑表达式 modulo 时pcop(m,r)减少到(-1)^m * r!, 所以,因为a = b + c(-1)取消,我们只剩下

choose(n,k) ≡ choose(a,b) * choose(A,B) (mod p)

在第二种情况下,我们发现

choose(n,k) = choose(a,b) * p * cop(a,A)/ (cop(b,B) * cop(c,C))

因为a = b + c + 1. 最后一位的进位意味着A < B,所以模p

p * cop(a,A) / (cop(b,B) * cop(c,C)) ≡ 0 = choose(A,B)

(我们可以用模逆的乘法代替除法,或者将其视为有理数的同余,这意味着分子可以被 整除p)。无论如何,我们再次发现

choose(n,k) ≡ choose(a,b) * choose(A,B) (mod p)

现在我们可以重现该choose(a,b)部分。

例子:

choose(144,6) (mod 5)
144 = 28 * 5 + 4
  6 =  1 * 5 + 1
choose(144,6) ≡ choose(28,1) * choose(4,1) (mod 5)
              ≡ choose(3,1) * choose(4,1) (mod 5)
              ≡ 3 * 4 = 12 ≡ 2 (mod 5)

choose(12349,789) ≡ choose(2469,157) * choose(4,4)
                  ≡ choose(493,31) * choose(4,2) * choose(4,4
                  ≡ choose(98,6) * choose(3,1) * choose(4,2) * choose(4,4)
                  ≡ choose(19,1) * choose(3,1) * choose(3,1) * choose(4,2) * choose(4,4)
                  ≡ 4 * 3 * 3 * 1 * 1 = 36 ≡ 1 (mod 5)

现在实现:

// Preconditions: 0 <= k <= n; p > 1 prime
long long choose_mod_one(long long n, long long k, long long p)
{
    // For small k, no recursion is necessary
    if (k < p) return choose_mod_two(n,k,p);
    long long q_n, r_n, q_k, r_k, choose;
    q_n = n / p;
    r_n = n % p;
    q_k = k / p;
    r_k = k % p;
    choose = choose_mod_two(r_n, r_k, p);
    // If the exponent of p in choose(n,k) isn't determined to be 0
    // before the calculation gets serious, short-cut here:
    /* if (choose == 0) return 0; */
    choose *= choose_mod_one(q_n, q_k, p);
    return choose % p;
}

// Preconditions: 0 <= k <= min(n,p-1); p > 1 prime
long long choose_mod_two(long long n, long long k, long long p)
{
    // reduce n modulo p
    n %= p;
    // Trivial checks
    if (n < k) return 0;
    if (k == 0 || k == n) return 1;
    // Now 0 < k < n, save a bit of work if k > n/2
    if (k > n/2) k = n-k;
    // calculate numerator and denominator modulo p
    long long num = n, den = 1;
    for(n = n-1; k > 1; --n, --k)
    {
        num = (num * n) % p;
        den = (den * k) % p;
    }
    // Invert denominator modulo p
    den = invert_mod(den,p);
    return (num * den) % p;
}

要计算模逆,可以使用费马(所谓的小)定理

如果p是素数a且不能被 整除p,则a^(p-1) ≡ 1 (mod p).

并将逆计算为a^(p-2) (mod p),或使用适用于更广泛参数范围的方法、扩展欧几里得算法或连分数展开式,它为您提供任何一对互质(正)整数的模逆:

long long invert_mod(long long k, long long m)
{
    if (m == 0) return (k == 1 || k == -1) ? k : 0;
    if (m < 0) m = -m;
    k %= m;
    if (k < 0) k += m;
    int neg = 1;
    long long p1 = 1, p2 = 0, k1 = k, m1 = m, q, r, temp;
    while(k1 > 0) {
        q = m1 / k1;
        r = m1 % k1;
        temp = q*p1 + p2;
        p2 = p1;
        p1 = temp;
        m1 = k1;
        k1 = r;
        neg = !neg;
    }
    return neg ? m - p2 : p2;
}

就像计算一样a^(p-2) (mod p),这是一种O(log p)算法,对于某些输入,它的速度要快得多(实际上是O(min(log k, log p)),因此对于 smallk和 large p,它要快得多),而对于其他输入则要慢得多。

总的来说,这种方式我们最多需要计算 O(log_p k) 个二项式系数 modulo p,其中每个二项式系数最多需要 O(p) 个操作,从而产生 O(p*log_p k) 个操作的总复杂度。当k显着大于时p,这比解要好得多O(k)。对于k <= p,它简化为O(k)具有一些开销的解决方案。

于 2012-06-02T13:59:20.873 回答
1

如果您不止一次计算它,还有另一种更快的方法。我将在 python 中发布代码,因为它可能是最容易转换成另一种语言的,尽管我会将 C++ 代码放在最后。

计算一次

蛮力:

def choose(n, k, m):
    ans = 1
    for i in range(k): ans *= (n-i)
    for i in range(k): ans //= i
    return ans % m

但是计算可能会涉及到非常大的数字,所以我们可以使用模块化的 airthmetic 技巧来代替:

(a * b) mod m = (a mod m) * (b mod m) mod m

(a / (b*c)) mod m = (a mod m) / ((b mod m) * (c mod m) mod m)

(a / b) mod m = (a mod m) * (b mod m)^-1

注意^-1最后一个等式末尾的 。这是bmod的乘法逆元m。它基本上意味着((b mod m) * (b mod m)^-1) mod m = 1,就像a * a^-1 = a * 1/a = 1使用(非零)整数一样。

这可以通过几种方式计算,其中之一是扩展欧几里得算法:

def multinv(n, m):
    ''' Multiplicative inverse of n mod m '''
    if m == 1: return 0
    m0, y, x = m, 0, 1

    while n > 1:
        y, x = x - n//m*y, y
        m, n = n%m, m
    
    return x+m0 if x < 0 else x

请注意,另一种方法,求幂,仅在m为素数时才有效。如果是,您可以这样做:

def powmod(b, e, m):
    ''' b^e mod m '''
    # Note: If you use python, there's a built-in pow(b, e, m) that's probably faster
    # But that's not in C++, so you can convert this instead:
    P = 1
    while e:
        if  e&1: P = P * b % m
        e >>= 1; b = b * b % m
    return P

def multinv(n, m):
    ''' Multiplicative inverse of n mod m, only if m is prime '''
    return powmod(n, m-2, m)
    

但请注意,扩展欧几里得算法往往仍然运行得更快,即使它们在技术上具有相同的时间复杂度 O(log m),因为它具有较低的常数因子。

所以现在完整的代码:

def multinv(n, m):
    ''' Multiplicative inverse of n mod m in log(m) '''
    if m == 1: return 0
    m0, y, x = m, 0, 1

    while n > 1:
        y, x = x - n//m*y, y
        m, n = n%m, m
    
    return x+m0 if x < 0 else x


def choose(n, k, m):
    num = den = 1
    for i in range(k): num = num * (n-i) % m
    for i in range(k): den = den * i % m
    return num * multinv(den, m)

多次查询

我们可以分别计算分子和分母,然后将它们组合起来。但请注意,我们为分子计算的乘积是n * (n-1) * (n-2) * (n-3) ... * (n-k+1)。如果您曾经了解过一种叫做前缀和的东西,那么这非常相似。所以让我们应用它。

预先fact[i] = i! mod m计算i最大值n,可能是1e7(一千万)。那么,分子是(fact[n] * fact[n-k]^-1) mod m,分母是fact[k]。所以我们可以计算choose(n, k, m) = fact[n] * multinv(fact[n-k], m) % m * multinv(fact[k], m) % m

Python代码:

MAXN = 1000 # Increase if necessary
MOD = 10**9+7 # A common mod that's used, change if necessary

fact = [1]
for i in range(1, MAXN+1):
    fact.append(fact[-1] * i % MOD)

def multinv(n, m):
    ''' Multiplicative inverse of n mod m in log(m) '''
    if m == 1: return 0
    m0, y, x = m, 0, 1

    while n > 1:
        y, x = x - n//m*y, y
        m, n = n%m, m
    
    return x+m0 if x < 0 else x


def choose(n, k, m):
    return fact[n] * multinv(fact[n-k] * fact[k] % m, m) % m

C++ 代码:

#include <iostream>
using namespace std;

const int MAXN = 1000; // Increase if necessary
const int MOD = 1e9+7; // A common mod that's used, change if necessary

int fact[MAXN+1];

int multinv(int n, int m) {
    /* Multiplicative inverse of n mod m in log(m) */
    if (m == 1) return 0;
    int m0 = m, y = 0, x = 1, t;

    while (n > 1) {
        t = y;
        y = x - n/m*y;
        x = t;
        
        t = m;
        m = n%m;
        n = t;
    }
    
    return x<0 ? x+m0 : x;
}

int choose(int n, int k, int m) {
    return (long long) fact[n]
         * multinv((long long) fact[n-k] * fact[k] % m, m) % m;
}

int main() {
    fact[0] = 1;
    for (int i = 1; i <= MAXN; i++) {
        fact[i] = (long long) fact[i-1] * i % MOD;
    }

    cout << choose(4, 2, MOD) << '\n';
    cout << choose(1e6, 1e3, MOD) << '\n';
}

请注意,我正在铸造long long以避免溢出。

于 2021-05-23T20:33:45.793 回答