C++ 计算大 n & k 的二项式系数 (nCk)

Calculating Binomial Coefficient (nCk) for large n & k

提问by

I just saw this question and have no idea how to solve it. can you please provide me with algorithms , C++ codes or ideas?

我刚看到这个问题，不知道如何解决。你能给我提供算法、C++ 代码或想法吗？

This is a very simple problem. Given the value of N and K, you need to tell us the value of the binomial coefficient C(N,K). You may rest assured that K <= N and the maximum value of N is 1,000,000,000,000,000. Since the value may be very large, you need to compute the result modulo 1009. Input

The first line of the input contains the number of test cases T, at most 1000. Each of the next T lines consists of two space separated integers N and K, where 0 <= K <= N and 1 <= N <= 1,000,000,000,000,000. Output

For each test case, print on a new line, the value of the binomial coefficient C(N,K) modulo 1009. Example

Input:
3
3 1
5 2
10 3

Output:
3
10
120

这是一个非常简单的问题。给定 N 和 K 的值，您需要告诉我们二项式系数 C(N,K) 的值。您可以放心，K <= N 并且 N 的最大值为 1,000,000,000,000,000。由于该值可能非常大，您需要计算结果模 1009。 输入

输入的第一行包含测试用例的数量 T，最多 1000。接下来的 T 行中的每一行都由两个空格分隔的整数 N 和 K 组成，其中 0 <= K <= N 和 1 <= N <= 1,000,000,000,000,000 . 输出

对于每个测试用例，在新行上打印二项式系数 C(N,K) 模 1009 的值。 示例

输入：
3
3 1
5 2
10 3

输出：
3
10
120

回答by

Notice that 1009 is a prime.

请注意，1009 是素数。

Now you can use Lucas' Theorem.

现在您可以使用卢卡斯定理。

Which states:

其中指出：

Let p be a prime.
If n  = a1a2...ar when written in base p and
if k  = b1b2...br when written in base p

(pad with zeroes if required)

Then

(n choose k) modulo p = (a1 choose b1) * (a2 choose  b2) * ... * (ar choose br) modulo p.

i.e. remainder of n choose k when divided by p is same as the remainder of
the product (a1 choose b1) * .... * (ar choose br) when divided by p.
Note: if bi > ai then ai choose bi is 0.

Thus your problem is reduced to finding the product modulo 1009 of at most log N/log 1009 numbers (number of digits of N in base 1009) of the form a choose b where a <= 1009 and b <= 1009.

因此，您的问题简化为找到最多 log N/log 1009 个数字（以 1009 为基数的 N 的位数）的乘积模 1009，其形式为 a 选择 b 其中 a <= 1009 和 b <= 1009。

This should make it easier even when N is close to 10^15.

即使 N 接近 10^15，这也应该使它更容易。

Note:

For N=10^15, N choose N/2 is more than 2^(100000000000000) which is way beyond an unsigned long long.

Also, the algorithm suggested by Lucas' theorem is O(log N) which is exponentiallyfaster than trying to compute the binomial coefficient directly (even if you did a mod 1009 to take care of the overflow issue).

笔记：

对于 N=10^15，N 选择 N/2 大于 2^(100000000000000)，这远远超出了 unsigned long long。

此外，卢卡斯定理建议的算法是 O(log N)，这 exponentially比尝试直接计算二项式系数要快（即使您使用 mod 1009 来处理溢出问题）。

Here is some code for Binomial I had written long back, all you need to do is to modify it to do the operations modulo 1009 (there might be bugs and not necessarily recommended coding style):

这是我很久以前写的二项式的一些代码，您需要做的就是修改它以执行取模 1009 的操作（可能存在错误，不一定是推荐的编码风格）：

class Binomial
{
public:
    Binomial(int Max)
    {
        max = Max+1;
        table = new unsigned int * [max]();
        for (int i=0; i < max; i++)
        {
            table[i] = new unsigned int[max]();

            for (int j = 0; j < max; j++)
            {
                table[i][j] = 0;
            }
        }
    }

    ~Binomial()
    {
        for (int i =0; i < max; i++)
        {
            delete table[i];
        }
        delete table;
    }

    unsigned int Choose(unsigned int n, unsigned int k);

private:
    bool Contains(unsigned int n, unsigned int k);

    int max;
    unsigned int **table;
};

unsigned int Binomial::Choose(unsigned int n, unsigned int k)
{
    if (n < k) return 0;
    if (k == 0 || n==1 ) return 1;
    if (n==2 && k==1) return 2;
    if (n==2 && k==2) return 1;
    if (n==k) return 1;


    if (Contains(n,k))
    {
        return table[n][k];
    }
    table[n][k] = Choose(n-1,k) + Choose(n-1,k-1);
    return table[n][k];
}

bool Binomial::Contains(unsigned int n, unsigned int k)
{
    if (table[n][k] == 0) 
    {
        return false;
    }
    return true;
}

回答by Steve Jessop

Binomial coefficient is one factorial divided by two others, although the k!term on the bottom cancels in an obvious way.

二项式系数是一个阶乘除以其他两个，尽管k!底部的项以明显的方式取消。

Observe that if 1009, (including multiples of it), appears more times in the numerator than the denominator, then the answer mod 1009 is 0. It can't appear more times in the denominator than the numerator (since binomial coefficients are integers), hence the only cases where you have to do anything are when it appears the same number of times in both. Don't forget to count multiples of (1009)^2 as two, and so on.

观察，如果 1009，（包括它的倍数）在分子中出现的次数多于分母，那么答案 mod 1009 是 0。它在分母中出现的次数不能多于分子（因为二项式系数是整数） ，因此您必须做任何事情的唯一情况是它在两者中出现的次数相同。不要忘记将 (1009)^2 的倍数计算为 2，依此类推。

After that, I think you're just mopping up small cases (meaning small numbers of values to multiply/divide), although I'm not sure without a few tests. On the plus side 1009 is prime, so arithmetic modulo 1009 takes place in a field, which means that after casting out multiples of 1009 from both top and bottom, you can do the rest of the multiplication and division mod 1009 in any order.

在那之后，我认为你只是在清理小案例（意味着要乘/除的少量值），尽管我不确定没有一些测试。从正面看，1009 是素数，因此算术模 1009 发生在一个域中，这意味着在从顶部和底部抛出 1009 的倍数后，您可以按任何顺序对其余的乘法和除法进行模 1009。

Where there are non-small cases left, they will still involve multiplying together long runs of consecutive integers. This can be simplified by knowing 1008! (mod 1009). It's -1 (1008 if you prefer), since 1 ... 1008 are the p-1non-zero elements of the prime field over p. Therefore they consist of 1, -1, and then (p-3)/2pairs of multiplicative inverses.

在剩余非小案例的情况下，它们仍将涉及将连续整数的长串相乘。这可以通过了解来简化1008! (mod 1009)。它是 -1（如果您愿意，则为 1008），因为 1 ... 1008 是p-1超过 的素数字段的非零元素p。因此，它们由 1、-1 和(p-3)/2成对的乘法逆组成。

So for example consider the case of C((1009^3), 200).

因此，例如考虑 C((1009^3), 200) 的情况。

Imagine that the number of 1009s are equal (don't know if they are, because I haven't coded a formula to find out), so that this is a case requiring work.

想象一下，1009的数量是相等的（不知道是不是，因为我没有写出公式来找出来），所以这是一个需要工作的案例。

On the top we have 201 ... 1008, which we'll have to calculate or look up in a precomputed table, then 1009, then 1010 ... 2017, 2018, 2019 ... 3026, 3027, etc. The ... ranges are all -1, so we just need to know how many such ranges there are.

在顶部，我们有 201 ... 1008，我们必须在预先计算的表中计算或查找，然后是 1009，然后是 1010 ... 2017、2018、2019 ... 3026、3027 等。 . .. 范围都是 -1，所以我们只需要知道有多少这样的范围。

That leaves 1009, 2018, 3027, which once we've cancelled them with 1009's from the bottom will just be 1, 2, 3, ... 1008, 1010, ..., plus some multiples of 1009^2, which again we'll cancel and leave ourselves with consecutive integers to multiply.

剩下 1009, 2018, 3027，一旦我们用底部的 1009 取消它们，它们将只是 1, 2, 3, ... 1008, 1010, ...，再加上一些 1009^2 的倍数，再次我们将取消并留下连续的整数进行乘法。

We can do something very similar with the bottom to compute the product mod 1009 of "1 ... 1009^3 - 200 with all the powers of 1009 divided out". That leaves us with a division in a prime field. IIRC that's tricky in principle, but 1009 is a small enough number that we can manage 1000 of them (the upper limit on the number of test cases).

我们可以用底部做一些非常相似的事情来计算“1 ... 1009^3 - 200 的乘积 mod 1009，其中 1009 的所有幂都被除掉了”。这给我们留下了一个主要领域的分裂。IIRC 原则上很棘手，但 1009 是一个足够小的数字，我们可以管理其中的 1000 个（测试用例数量的上限）。

Of course with k=200, there's an enormous overlap which could be cancelled more directly. That's what I meant by small cases and non-small cases: I've treated it like a non-small case, when in fact we could get away with just "brute-forcing" this one, by calculating ((1009^3-199) * ... * 1009^3) / 200!

当然，在 k=200 的情况下，存在可以更直接取消的巨大重叠。这就是我所说的小案例和非小案例的意思：我把它当作一个非小案例来对待，而实际上我们可以通过计算“蛮力”来摆脱这个问题((1009^3-199) * ... * 1009^3) / 200!

回答by dmuir

I don't think you want to calculate C(n,k) and then reduce mod 1009. The biggest one, C(1e15,5e14) will require something like 1e16 bits ~ 1000 terabytes

我不认为你想计算 C(n,k) 然后减少 mod 1009。最大的，C(1e15,5e14) 将需要类似 1e16 位 ~ 1000 TB 的东西

Moreover executing the loop in snakiles answer 1e15 times seems like it might take a while. What you might use is, if

此外，在 snakiles 中执行循环 1e15 次似乎可能需要一段时间。你可能会使用的是，如果

n = n0 + n1*p + n2*p^2 ... + nd*p^d

n = n0 + n1*p + n2*p^2 ... + nd*p^d

m = m0 + m1*p + m2*p^2 ... + md*p^d

m = m0 + m1*p + m2*p^2 ... + md*p^d

(where 0<=mi,ni < p)

（其中 0<=mi,ni < p）

then C(n,m) = C(n0,m0) * C(n1,m1) *... * C(nd, nd) mod p

然后 C(n,m) = C(n0,m0) * C(n1,m1) *... * C(nd, nd) mod p

see, eg http://www.cecm.sfu.ca/organics/papers/granville/paper/binomial/html/binomial.html

参见，例如http://www.cecm.sfu.ca/organics/papers/granville/paper/binomial/html/binomial.html

One way would be to use pascal's triangle to build a table of all C(m,n) for 0<=m<=n<=1009.

一种方法是使用帕斯卡三角形为 0<=m<=n<=1009 构建一个包含所有 C(m,n) 的表。

回答by snakile

psudo code for calculating nCk:

用于计算 nCk 的伪代码：

result = 1    
for i=1 to min{K,N-K}:
   result *= N-i+1
   result /= i
return result

Time Complexity: O(min{K,N-K})

时间复杂度：O(min{K,NK})

The loop goes from i=1 to min{K,N-K}instead of from i=1 to K, and that's ok because

循环from i=1 to min{K,N-K}代替from i=1 to K，这没关系，因为

C(k,n) = C(k, n-k)

And you can calculate the thing even more efficiently if you use the GammaLn function.

如果您使用 GammaLn 函数，您可以更有效地计算事物。

nCk = exp(GammaLn(n+1)-GammaLn(k+1)-GammaLn(n-k+1))

The GammaLn function is the natural logarithm of the Gamma function. I know there's an efficient algorithm to calculate the GammaLn function but that algorithm isn't trivial at all.

GammaLn 函数是Gamma 函数的自然对数。我知道有一种有效的算法来计算 GammaLn 函数，但该算法一点也不简单。

回答by AndyUpNorth

The following code shows how to obtain all the binomial coefficients for a given size 'n'. You could easily modify it to stop at a given k in order to determine nCk. It is computationally very efficient, it's simple to code, and works for very large n and k.

以下代码显示了如何获取给定大小“n”的所有二项式系数。您可以轻松修改它以在给定的 k 处停止，以确定 nCk。它在计算上非常高效，编码简单，并且适用于非常大的 n 和 k。

binomial_coefficient = 1
output(binomial_coefficient)
col = 0
n = 5

do while col < n
    binomial_coefficient = binomial_coefficient * (n + 1 - (col + 1)) / (col + 1)
    output(binomial_coefficient)
    col = col + 1
loop

The output of binomial coefficients is therefore:

因此二项式系数的输出是：

1
1 *  (5 + 1 - (0 + 1)) / (0 + 1) = 5 
5 *  (5 + 1 - (1 + 1)) / (1 + 1) = 15
15 * (5 + 1 - (2 + 1)) / (2 + 1) = 15 
15 * (5 + 1 - (3 + 1)) / (3 + 1) = 5 
5 *  (5 + 1 - (4 + 1)) / (4 + 1) = 1

I had found the formula once upon a time on Wikipedia but for some reason it's no longer there :(

我曾在维基百科上找到过这个公式，但由于某种原因它不再存在:(