n选k mod p ?

时间：2021-02-19 15:01:56

What I mean by "large n" is something in the millions. p is prime.

我所说的“大n”指的是数以百万计的人。p是质数。

I've tried http://apps.topcoder.com/wiki/display/tc/SRM+467 But the function seems to be incorrect (I tested it with 144 choose 6 mod 5 and it gives me 0 when it should give me 2)

我已经尝试了http://apps.topcoder.com/wiki/display/tc/SRM+467，但是这个函数似乎是不正确的(我用144选择了6 mod 5，它给了我0，它应该给我2)

I've tried http://online-judge.uva.es/board/viewtopic.php?f=22&t=42690 But I don't understand it fully

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

I've also made a memoized recursive function that uses the logic (combinations(n-1, k-1, p)%p + combinations(n-1, k, p)%p) but it gives me stack overflow problems because n is large

我还做了一个memoized递归函数，它使用了逻辑(组合(n-1, k-1, p)%p +组合(n-1, k, p)%p)，但它给了我堆栈溢出问题，因为n很大。

I've tried Lucas Theorem but it appears to be either slow or inaccurate.

我已经试过了卢卡斯定理，但它看起来要么慢要么不准确。

All I'm trying to do is create a fast/accurate n choose k mod p for large n. If anyone could help show me a good implementation for this I'd be very grateful. Thanks.

我所要做的就是创建一个快速/准确的n选k对大n，如果有人能帮我展示一个好的实现，我将非常感激。谢谢。

As requested, the memoized version that hits stack overflows for large n:

按要求，对大n的堆栈溢出的memoized版本:

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;
   }  
}

2 个解决方案

#1

51  

So, here is how you can solve your problem.

所以，这里是你如何解决你的问题。

Of course you know the formula:

你当然知道公式:

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

(See http://en.wikipedia.org/wiki/Binomial_coefficient#Computing_the_value_of_binomial_coefficients)

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

You know how to compute the numerator:

你知道如何计算分子:

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

Now, as p is prime the reciprocal of each integer that is coprime with p is well defined i.e. a-1 can be found. And this can be done using Fermat's theorem ap-1=1(mod p) => a*ap-2=1(mod p) and so a-1=ap-2. Now all you need to do is to implement fast exponentiation(for example using the binary method):

现在，当p是质数时，每一个整数的倒数都被定义为a-1。这可以用费马定理ap-1=1(mod p) => a*ap-2=1(mod p)和a-1=ap-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;
}

And now you can add the denominator to our result:

现在你可以把分母加入到我们的结果中:

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

Please note I am using long long everywhere to avoid type overflow. Of course you don't need to do k exponentiations - you can compute k!(mod p) and then divide only once:

请注意，我使用了很长时间，以避免类型溢出。当然你不需要做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;

EDIT: as per @dbaupp's comment if k >= p the k! will be equal to 0 modulo p and (k!)^-1 will not be defined. To avoid that first compute the degree with which p is in n*(n-1)...(n-k+1) and in k! and compare them:

编辑:根据@dbaupp的评论，如果k >= p 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;
}

EDIT: There is one more optimization that can be added to the solution above - instead of computing the inverse number of each multiple in k!, we can compute k!(mod p) and then compute the inverse of that number. Thus we have to pay the logarithm for the exponentiation only once. Of course again we have to discard the p divisors of each multiple. We only have to change the last loop with this:

编辑:还有一个优化可以添加到上面的解决方案中——而不是计算k中的每个倍数的逆数!我们可以计算k!(mod p)然后计算这个数的倒数。因此，我们只需要支付一次取幂的对数。当然，我们还必须抛弃每个倍数的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;

#2

12  

For large k, we can reduce the work significantly by exploiting two fundamental facts:

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

  1. If p is a prime, the exponent of p in the prime factorisation of n! is given by (n - s_p(n)) / (p-1), where s_p(n) is the sum of the digits of n in the base p representation (so for p = 2, it's popcount). Thus the exponent of p in the prime factorisation of choose(n,k) is (s_p(k) + s_p(n-k) - s_p(n)) / (p-1), in particular, it is zero if and only if the addition k + (n-k) has no carry when performed in base p (the exponent is the number of carries).

    如果p是质数，那么p在n的主要因子中的指数。是(n - s_p(n)) / (p-1)，其中s_p(n)是n在基本p表示中的数字之和(因此p = 2，它是popcount)。因此，在选择(n,k)的主要因子中p的指数为(s_p(k) + s_p(n-k) - s_p(n)) / (p-1)，特别是当且仅当k + (n-k)在基底p中执行时没有进位(指数为进位数)。

  2. Wilson's theorem: p is a prime, if and only if (p-1)! ≡ (-1) (mod p).

    威尔逊定理:p是一个质数，如果且仅当(p-1)!≡(1)(mod p)。

The exponent of p in the factorisation of n! is usually calculated by

在n的因子化中p的指数!通常是计算

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

The check for divisibility of choose(n,k) by p is not strictly necessary, but it's reasonable to have that first, since it will often be the case, and then it's less work:

p (n,k)的可分度检查并不是严格意义上的，但它是合理的，因为它经常是这样的，然后它的工作就更少了:

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

Now let us take a closer look at n!. We separate the numbers ≤ n into the multiples of p and the numbers coprime to p. With

现在让我们仔细看看n!我们把这些数分离成p的倍数和coprime到p的数。

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

The multiples of p contribute p^q * q!. The numbers coprime to p contribute the product of (j*p + k), 1 ≤ k < p for 0 ≤ j < q, and the product of (q*p + k), 1 ≤ k ≤ r.

的倍数p p ^ *贡献!。数字coprime p贡献的产物(j * p + k),1≤0≤j k < p < q(q * p + k)的产品,1 k≤≤r。

For the numbers coprime to p we will only be interested in the contribution modulo p. Each of the full runs j*p + k, 1 ≤ k < p is congruent to (p-1)! modulo p, so altogether they produce a contribution of (-1)^q modulo p. The last (possibly) incomplete run produces r! modulo p.

数字coprime p我们才会感兴趣的贡献模p。每一个完整的运行j * p + k,1≤k < p是一致的(p - 1)!模p，所以它们总共产生了(-1)q模p。最后(可能)不完整的运行产生r!模p。

So if we write

如果我们写

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

we get

我们得到了

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

where cop(m,r) is the product of all numbers coprime to p which are ≤ m*p + r.

警察(m,r)是所有数字的乘积coprime p≤m * p + r。

There are two possibilities, a = b + c and A = B + C, or a = b + c + 1 and A = B + C - p.

有两种可能，a = b + c, a = b + c, a = b + c + 1, a = b + c - p。

In our calculation, we have eliminated the second possibility beforehand, but that is not essential.

在我们的计算中，我们已经提前排除了第二种可能性，但这并不重要。

In the first case, the explicit powers of p cancel, and we are left with

在第一种情况下，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))

Any powers of p dividing choose(n,k) come from choose(a,b) - in our case, there will be none, since we've eliminated these cases before - and, although cop(a,A) / (cop(b,B) * cop(c,C)) need not be an integer (consider e.g. choose(19,9) (mod 5)), when considering the expression modulo p, cop(m,r) reduces to (-1)^m * r!, so, since a = b + c, the (-1) cancel and we are left with

任何权力的p分选择(n,k)来自选择(a,b)——在我们的例子中,将会有,因为我们之前消除这些情况下,尽管警察(a)/(cop(b,b)*警察(c,c))不需要一个整数(考虑如选择(19日9)5(mod)),在考虑表达模p,警察(m,r)减少(1)r m ^ * !因为a = b + c，(-1)消掉了，剩下了。

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

In the second case, we find

在第二个例子中，我们发现。

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

since a = b + c + 1. The carry in the last digit means that A < B, so modulo p

因为a = b + c + 1。最后一个数字的进位意味着A < B，所以模p。

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

(where we can either replace the division with a multiplication by the modular inverse, or view it as a congruence of rational numbers, meaning the numerator is divisible by p). Anyway, we again find

(在这里，我们可以用模块化的逆矩阵来替换这个除法，或者把它看成是一个合理的数字的同余，这意味着分子可以被p整除)。无论如何，我们再次找到了。

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

Now we can recur for the choose(a,b) part.

现在我们可以重新选择(a,b)部分。

Example:

例子:

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)

Now the implementation:

现在的实现:

// 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;
}

To calculate the modular inverse, you can use Fermat's (so-called little) theorem

为了计算模块的逆，你可以使用费马的(所谓的小)定理。

If p is prime and a not divisible by p, then a^(p-1) ≡ 1 (mod p).

如果p '和一个不是p整除,那么^(p - 1)≡1 p(mod)。

and calculate the inverse as a^(p-2) (mod p), or use a method applicable to a wider range of arguments, the extended Euclidean algorithm or continued fraction expansion, which give you the modular inverse for any pair of coprime (positive) integers:

计算逆矩阵(p-2) (p-2) (mod p)，或者使用一种适用于更广范围的参数的方法，扩展的欧几里得算法或继续分式展开，它给你的任意一对coprime(正整数)的模逆:

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;
}

Like calculating a^(p-2) (mod p), this is an O(log p) algorithm, for some inputs it's significantly faster (it's actually O(min(log k, log p)), so for small k and large p, it's considerably faster), for others it's slower.

像计算^(p 2)(mod p),这是一个O(log p)算法,对一些输入的更快(实际上是O(min(对数k,p)),所以对于小k和大p,它相当快),为别人慢。

Overall, this way we need to calculate at most O(log_p k) binomial coefficients modulo p, where each binomial coefficient needs at most O(p) operations, yielding a total complexity of O(p*log_p k) operations. When k is significantly larger than p, that is much better than the O(k) solution. For k <= p, it reduces to the O(k) solution with some overhead.

总而言之，我们需要计算最多O(log_p k)二项式系数，其中每个二项式系数在大多数O(p)操作中都需要，从而产生O(p*log_p k)操作的总复杂度。当k明显大于p时，这比O(k)解要好得多。对于k <= p，它通过一些开销减少到O(k)解决方案。

#1

51  

So, here is how you can solve your problem.

所以，这里是你如何解决你的问题。

Of course you know the formula:

你当然知道公式:

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

(See http://en.wikipedia.org/wiki/Binomial_coefficient#Computing_the_value_of_binomial_coefficients)

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

You know how to compute the numerator:

你知道如何计算分子:

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

Now, as p is prime the reciprocal of each integer that is coprime with p is well defined i.e. a-1 can be found. And this can be done using Fermat's theorem ap-1=1(mod p) => a*ap-2=1(mod p) and so a-1=ap-2. Now all you need to do is to implement fast exponentiation(for example using the binary method):

现在，当p是质数时，每一个整数的倒数都被定义为a-1。这可以用费马定理ap-1=1(mod p) => a*ap-2=1(mod p)和a-1=ap-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;
}

And now you can add the denominator to our result:

现在你可以把分母加入到我们的结果中:

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

Please note I am using long long everywhere to avoid type overflow. Of course you don't need to do k exponentiations - you can compute k!(mod p) and then divide only once:

请注意，我使用了很长时间，以避免类型溢出。当然你不需要做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;

EDIT: as per @dbaupp's comment if k >= p the k! will be equal to 0 modulo p and (k!)^-1 will not be defined. To avoid that first compute the degree with which p is in n*(n-1)...(n-k+1) and in k! and compare them:

编辑:根据@dbaupp的评论，如果k >= p 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;
}

EDIT: There is one more optimization that can be added to the solution above - instead of computing the inverse number of each multiple in k!, we can compute k!(mod p) and then compute the inverse of that number. Thus we have to pay the logarithm for the exponentiation only once. Of course again we have to discard the p divisors of each multiple. We only have to change the last loop with this:

编辑:还有一个优化可以添加到上面的解决方案中——而不是计算k中的每个倍数的逆数!我们可以计算k!(mod p)然后计算这个数的倒数。因此，我们只需要支付一次取幂的对数。当然，我们还必须抛弃每个倍数的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;

#2

12  

For large k, we can reduce the work significantly by exploiting two fundamental facts:

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

  1. If p is a prime, the exponent of p in the prime factorisation of n! is given by (n - s_p(n)) / (p-1), where s_p(n) is the sum of the digits of n in the base p representation (so for p = 2, it's popcount). Thus the exponent of p in the prime factorisation of choose(n,k) is (s_p(k) + s_p(n-k) - s_p(n)) / (p-1), in particular, it is zero if and only if the addition k + (n-k) has no carry when performed in base p (the exponent is the number of carries).

    如果p是质数，那么p在n的主要因子中的指数。是(n - s_p(n)) / (p-1)，其中s_p(n)是n在基本p表示中的数字之和(因此p = 2，它是popcount)。因此，在选择(n,k)的主要因子中p的指数为(s_p(k) + s_p(n-k) - s_p(n)) / (p-1)，特别是当且仅当k + (n-k)在基底p中执行时没有进位(指数为进位数)。

  2. Wilson's theorem: p is a prime, if and only if (p-1)! ≡ (-1) (mod p).

    威尔逊定理:p是一个质数，如果且仅当(p-1)!≡(1)(mod p)。

The exponent of p in the factorisation of n! is usually calculated by

在n的因子化中p的指数!通常是计算

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

The check for divisibility of choose(n,k) by p is not strictly necessary, but it's reasonable to have that first, since it will often be the case, and then it's less work:

p (n,k)的可分度检查并不是严格意义上的，但它是合理的，因为它经常是这样的，然后它的工作就更少了:

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

Now let us take a closer look at n!. We separate the numbers ≤ n into the multiples of p and the numbers coprime to p. With

现在让我们仔细看看n!我们把这些数分离成p的倍数和coprime到p的数。

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

The multiples of p contribute p^q * q!. The numbers coprime to p contribute the product of (j*p + k), 1 ≤ k < p for 0 ≤ j < q, and the product of (q*p + k), 1 ≤ k ≤ r.

的倍数p p ^ *贡献!。数字coprime p贡献的产物(j * p + k),1≤0≤j k < p < q(q * p + k)的产品,1 k≤≤r。

For the numbers coprime to p we will only be interested in the contribution modulo p. Each of the full runs j*p + k, 1 ≤ k < p is congruent to (p-1)! modulo p, so altogether they produce a contribution of (-1)^q modulo p. The last (possibly) incomplete run produces r! modulo p.

数字coprime p我们才会感兴趣的贡献模p。每一个完整的运行j * p + k,1≤k < p是一致的(p - 1)!模p，所以它们总共产生了(-1)q模p。最后(可能)不完整的运行产生r!模p。

So if we write

如果我们写

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

we get

我们得到了

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

where cop(m,r) is the product of all numbers coprime to p which are ≤ m*p + r.

警察(m,r)是所有数字的乘积coprime p≤m * p + r。

There are two possibilities, a = b + c and A = B + C, or a = b + c + 1 and A = B + C - p.

有两种可能，a = b + c, a = b + c, a = b + c + 1, a = b + c - p。

In our calculation, we have eliminated the second possibility beforehand, but that is not essential.

在我们的计算中，我们已经提前排除了第二种可能性，但这并不重要。

In the first case, the explicit powers of p cancel, and we are left with

在第一种情况下，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))

Any powers of p dividing choose(n,k) come from choose(a,b) - in our case, there will be none, since we've eliminated these cases before - and, although cop(a,A) / (cop(b,B) * cop(c,C)) need not be an integer (consider e.g. choose(19,9) (mod 5)), when considering the expression modulo p, cop(m,r) reduces to (-1)^m * r!, so, since a = b + c, the (-1) cancel and we are left with

任何权力的p分选择(n,k)来自选择(a,b)——在我们的例子中,将会有,因为我们之前消除这些情况下,尽管警察(a)/(cop(b,b)*警察(c,c))不需要一个整数(考虑如选择(19日9)5(mod)),在考虑表达模p,警察(m,r)减少(1)r m ^ * !因为a = b + c，(-1)消掉了，剩下了。

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

In the second case, we find

在第二个例子中，我们发现。

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

since a = b + c + 1. The carry in the last digit means that A < B, so modulo p

因为a = b + c + 1。最后一个数字的进位意味着A < B，所以模p。

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

(where we can either replace the division with a multiplication by the modular inverse, or view it as a congruence of rational numbers, meaning the numerator is divisible by p). Anyway, we again find

(在这里，我们可以用模块化的逆矩阵来替换这个除法，或者把它看成是一个合理的数字的同余，这意味着分子可以被p整除)。无论如何，我们再次找到了。

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

Now we can recur for the choose(a,b) part.

现在我们可以重新选择(a,b)部分。

Example:

例子:

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)

Now the implementation:

现在的实现:

// 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;
}

To calculate the modular inverse, you can use Fermat's (so-called little) theorem

为了计算模块的逆，你可以使用费马的(所谓的小)定理。

If p is prime and a not divisible by p, then a^(p-1) ≡ 1 (mod p).

如果p '和一个不是p整除,那么^(p - 1)≡1 p(mod)。

and calculate the inverse as a^(p-2) (mod p), or use a method applicable to a wider range of arguments, the extended Euclidean algorithm or continued fraction expansion, which give you the modular inverse for any pair of coprime (positive) integers:

计算逆矩阵(p-2) (p-2) (mod p)，或者使用一种适用于更广范围的参数的方法，扩展的欧几里得算法或继续分式展开，它给你的任意一对coprime(正整数)的模逆:

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;
}

Like calculating a^(p-2) (mod p), this is an O(log p) algorithm, for some inputs it's significantly faster (it's actually O(min(log k, log p)), so for small k and large p, it's considerably faster), for others it's slower.

像计算^(p 2)(mod p),这是一个O(log p)算法,对一些输入的更快(实际上是O(min(对数k,p)),所以对于小k和大p,它相当快),为别人慢。

Overall, this way we need to calculate at most O(log_p k) binomial coefficients modulo p, where each binomial coefficient needs at most O(p) operations, yielding a total complexity of O(p*log_p k) operations. When k is significantly larger than p, that is much better than the O(k) solution. For k <= p, it reduces to the O(k) solution with some overhead.

总而言之，我们需要计算最多O(log_p k)二项式系数，其中每个二项式系数在大多数O(p)操作中都需要，从而产生O(p*log_p k)操作的总复杂度。当k明显大于p时，这比O(k)解要好得多。对于k <= p，它通过一些开销减少到O(k)解决方案。

标签：

相关文章