logo标签列表

计算组合数量

c++algorithmcombinatorics
31

您好，

我知道您可以使用以下公式（不考虑重复且顺序不重要）获取组合数：

//从n个物品中选r个
n! / r!(n - r)!

但是，我不知道如何在C++中实现它，因为例如：

n = 52
n! = 8,0658175170943878571660636856404e+67

即使对于 unsigned __int64 或 unsigned long long 类型的变量而言，这个数字也太大了。请问有没有不需要第三方“bigint”库的解决方法？

有关使用大整数，请参见以下链接：<br> https://dev59.com/5nM_5IYBdhLWcg3w3nXz<br> https://dev59.com/VXNA5IYBdhLWcg3wL62x<br> https://dev59.com/cnVC5IYBdhLWcg3wnCf6<br> - Sajad Bahmani
11个回答
43

这是一个古老的算法，精确无误，除非结果太大而无法容纳在long long中，否则不会溢出。

unsigned long long
choose(unsigned long long n, unsigned long long k) {
    if (k > n) {
        return 0;
    }
    unsigned long long r = 1;
    for (unsigned long long d = 1; d <= k; ++d) {
        r *= n--;
        r /= d;
    }
    return r;
}

这个算法也可以在Knuth的《计算机程序设计艺术》第3版第2卷：数字半数学算法中找到。

更新：有一种小可能性，即算法会在以下行上溢出：

r *= n--;

对于非常大的n，一个天真的上界是sqrt(std::numeric_limits<long long>::max())，这意味着n小于大约4,000,000,000。

这个可以通过 r *= (n--) / d 来改进，先进行除法运算吗？ - GManNickG
2
GManNickG，在我看来，这样做会失去精度。 - Andres Riofrio
2
一种改进方法是将k设置为k和（n-k）的最小值。 - Agrim Pathak
2
正如Howard在这里所示，这个答案不够精确，特别是从for very large n. A naive upper bound...开始。 - Antonio
那么“update”可以说是非常不正确的吗？（下一个被投票的答案显示n == 67时发生了溢出）。 - Will Ness
35

以下内容来自Andreas的回答:

Here's an ancient algorithm which is exact and doesn't overflow unless the result is to big for a long long

unsigned long long
choose(unsigned long long n, unsigned long long k) {
    if (k > n) {
        return 0;
    }
    unsigned long long r = 1;
    for (unsigned long long d = 1; d <= k; ++d) {
        r *= n--;
        r /= d;
    }
    return r;
}

This algorithm is also in Knuth's "The Art of Computer Programming, 3rd Edition, Volume 2: Seminumerical Algorithms" I think.

UPDATE: There's a small possibility that the algorithm will overflow on the line:

r *= n--;

for very large n. A naive upper bound is sqrt(std::numeric_limits<long long>::max()) which means an n less than rougly 4,000,000,000.

考虑 n == 67 和 k == 33 的情况。上述算法在使用64位无符号长整型时会溢出。然而正确答案可以用64位表示：14,226,520,737,620,288,370。上述算法对其溢出保持沉默，choose(67, 33) 返回：
8,829,174,638,479,413
这是一个看起来可信但不正确的答案。
然而，只要最终答案可表示，上述算法可以稍加修改而不会溢出。
关键在于认识到每次迭代时，除法 r/d 是准确的。暂时改写为：
r = r * n / d;
--n;

准确来说，如果您将r、n和d扩展为它们的质因数分解，那么可以轻松取消d，并留下一个修改后的值n，称其为t，然后计算r就是：

// compute t from r, n and d
r = r * t;
--n;

一种快速简便的方法是找到r和d的最大公约数，称其为g:

unsigned long long g = gcd(r, d);
// now one can divide both r and d by g without truncation
r /= g;
unsigned long long d_temp = d / g;
--n;

现在我们可以使用d_temp和n（找到最大公约数）做同样的事情。但是，由于我们事先知道r * n / d是精确的，因此我们也知道gcd（d_temp，n）== d_temp，因此我们不需要计算它。因此，我们可以通过d_temp将n除以：
unsigned long long g = gcd(r, d);
// now one can divide both r and d by g without truncation
r /= g;
unsigned long long d_temp = d / g;
// now one can divide n by d/g without truncation
unsigned long long t = n / d_temp;
r = r * t;
--n;

清理：

unsigned long long
gcd(unsigned long long x, unsigned long long y)
{
    while (y != 0)
    {
        unsigned long long t = x % y;
        x = y;
        y = t;
    }
    return x;
}

unsigned long long
choose(unsigned long long n, unsigned long long k)
{
    if (k > n)
        throw std::invalid_argument("invalid argument in choose");
    unsigned long long r = 1;
    for (unsigned long long d = 1; d <= k; ++d, --n)
    {
        unsigned long long g = gcd(r, d);
        r /= g;
        unsigned long long t = n / (d / g);
        if (r > std::numeric_limits<unsigned long long>::max() / t)
           throw std::overflow_error("overflow in choose");
        r *= t;
    }
    return r;
}

现在您可以计算choose（67，33）而不会溢出。如果尝试计算choose（68，33），您将收到异常而不是错误答案。

Howard，我已经修复了你回答中混乱的格式。请阅读编辑窗格右侧的编辑提示，以了解如何自行完成此操作。哦，欢迎来到SO！ - sbi
1
@sbi想引用被接受的答案，这就是为什么看起来有点奇怪。公正地说，编辑对某些事情真的很糟糕。 - Johannes Schaub - litb
@Johannes：哦，我完全没注意到！也许提示会更合适？ - sbi
你的编辑非常准确，非常感谢！我在这里还是个新手，正在学习适当的礼仪和编辑技巧。 - Howard Hinnant
同样的优化方案也可以应用于这里：将k设置为k和n-k中的最小值... - Aconcagua
6
下面的例程将使用递归定义和备忘录计算n选k。该例程非常快速和准确：
inline unsigned long long n_choose_k(const unsigned long long& n,
                                     const unsigned long long& k)
{
   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;       
   typedef unsigned long long value_type;
   value_type* table = new value_type[static_cast<std::size_t>(n * n)];
   std::fill_n(table,n * n,0);
   class n_choose_k_impl
   {
   public:

      n_choose_k_impl(value_type* table,const value_type& dimension)
      : table_(table),
        dimension_(dimension)
      {}

      inline value_type& lookup(const value_type& n, const value_type& k)
      {
         return table_[dimension_ * n + k];
      }

      inline value_type compute(const value_type& n, const value_type& k)
      {
         if ((0 == k) || (k == n))
            return 1;
         value_type v1 = lookup(n - 1,k - 1);
         if (0 == v1)
            v1 = lookup(n - 1,k - 1) = compute(n - 1,k - 1);
         value_type v2 = lookup(n - 1,k);
         if (0 == v2)
            v2 = lookup(n - 1,k) = compute(n - 1,k);
         return v1 + v2;
      }

      value_type* table_;
      value_type dimension_;
   };
   value_type result = n_choose_k_impl(table,n).compute(n,k);
   delete [] table;
   return result;
}
4
请记住：

n! / ( n - r )! = n * ( n - 1) * .. * (n - r + 1 )

因此，它比n！要小得多。所以解决方案是计算n *（n-1）* ... *（n-r + 1），而不是先计算n！然后再除以它。
当然，这完全取决于n和r的相对大小-如果r相对于n很大，则仍然无法适应。
请注意，问题是如何计算n! / r!(n - r)!，而不是n! / (n - r)!。 - wenqiang
在你的回答中，好像缺少了除以 r! 的部分，你似乎只计算了 n!/(n-r)!。 - Antonio
2

好的，我必须回答自己的问题。我正在阅读关于帕斯卡三角形的内容，并偶然发现我们可以用它来计算组合的数量：

#include <iostream>
#include <boost/cstdint.hpp>

boost::uint64_t Combinations(unsigned int n, unsigned int r)
{
    if (r > n)
        return 0;

    /** We can use Pascal's triange to determine the amount
      * of combinations. To calculate a single line:
      *
      * v(r) = (n - r) / r
      *
      * Since the triangle is symmetrical, we only need to calculate
      * until r -column.
      */

    boost::uint64_t v = n--;

    for (unsigned int i = 2; i < r + 1; ++i, --n)
        v = v * n / i;

    return v;
}

int main()
{
    std::cout << Combinations(52, 5) << std::endl;
}
1
是的，这正是我发布的相同算法。恭喜你自己想出来了 ;) - Andreas Brinck
注意：自从C++11，uint64_t已经成为#include <cstdint>的一部分，因此我们不再需要在此示例中使用boost。 - M.M
1
获取二项式系数的质因数分解可能是计算它最有效的方法，特别是如果乘法很昂贵。这在计算阶乘的相关问题中肯定是正确的（例如，请参见此处）。这里有一个基于埃拉托斯特尼筛法的简单算法来计算质因数分解。基本思想是通过使用筛法找到质数，但同时也计算它们的倍数中有多少个落在[1，k]和[n-k+1，n]范围内。筛法本质上是一个O(n log log n)算法，但没有执行乘法。一旦找到质因数分解，实际所需的乘法次数最坏为O（n log log n / log n），可能还有更快的方法。
prime_factors = []

n = 20
k = 10

composite = [True] * 2 + [False] * n

for p in xrange(n + 1):
if composite[p]:
    continue

q = p
m = 1
total_prime_power = 0
prime_power = [0] * (n + 1)

while True:

    prime_power[q] = prime_power[m] + 1
    r = q

    if q <= k:
        total_prime_power -= prime_power[q]

    if q > n - k:
        total_prime_power += prime_power[q]

    m += 1
    q += p

    if q > n:
        break

    composite[q] = True

prime_factors.append([p, total_prime_power])

 print prime_factors
1

对 Howard Hinnant 在这个问题中的回答进行了一点改进： 每次循环调用 gcd() 似乎有点慢。 我们可以将 gcd() 调用聚合到最后一个调用中，同时充分利用 Knuth 的书《计算机程序设计艺术》第三版第二卷中的标准算法：

const uint64_t u64max = std::numeric_limits<uint64_t>::max();
uint64_t choose(uint64_t n, uint64_t k)
{
    if (k > n)
        throw std::invalid_argument(std::string("invalid argument in ") + __func__);

    if (k > n - k)
        k = n - k;

    uint64_t r = 1;
    uint64_t d;
    for (d = 1; d <= k; ++d) {
        if (r > u64max / n)
            break;
        r *= n--;
        r /= d;
    }

    if (d > k)
        return r;

    // Let N be the original n,
    // n is the current n (when we reach here)
    // We want to calculate C(N,k),
    // Currently we already calculated the r value so far:
    // r = C(N, n) = C(N, N-n) = C(N, d-1)
    // Note that N-n = d-1
    // In addition we know the following identity formula:
    //  C(N,k) = C(N,d-1) * C(N-d+1, k-d+1) / C(k, k-d+1)
    //         = C(N,d-1) * C(n, k-d+1) / C(k, k-d+1)
    // Using this formula, we effectively reduce the calculation,
    // while recursively use the same function.
    uint64_t b = choose(n, k-d+1);
    if (b == u64max) {
        return u64max;  // overflow
    }

    uint64_t c = choose(k, k-d+1);
    if (c == u64max) {
        return u64max;  // overflow
    }

    // Now, the combinatorial should be r * b / c
    // We can use gcd() to calculate this:
    // We Pick b for gcd: b < r almost (if not always) in all cases
    uint64_t g = gcd(b, c);
    b /= g;
    c /= g;
    r /= c;

    if (r > u64max / b)
        return u64max;   // overflow

    return r * b;
}

请注意，递归深度通常为2（我真的没有看到需要3次调用choose()的情况，组合减少得相当不错），即对于非溢出情况调用choose() 3次。
如果您喜欢，可以使用unsigned long long替换uint64_t。
也许更快的替代方法是使用大数乘法/除法来计算“r x b / c”的最后一部分。 - Robin Hsu
1

使用一个长双精度浮点数的恶意技巧，可以获得与Howard Hinnant相同甚至更高的精度：

unsigned long long n_choose_k(int n, int k)
{
    long double f = n;
    for (int i = 1; i<k+1; i++)
        f /= i;
    for (int i=1; i<k; i++)
        f *= n - i;

    unsigned long long f_2 = std::round(f);

    return f_2;
}

这个想法是先除以k!，然后再乘以n(n-1)...(n-k+1)。通过改变for循环的顺序，可以避免使用double进行近似计算。
0

最短的方法之一：

int nChoosek(int n, int k){
    if (k > n) return 0;
    if (k == 0) return 1;
    return nChoosek(n - 1, k) + nChoosek(n - 1, k - 1);
}
0
一个类似于埃拉托斯特尼筛法的方法。 虽然埃拉托斯特尼筛法是一种多重消除的方法， 但这个方法是一种多重半消除的方法。 由于n！/((n-r)!r!)始终为整数，首先取消分母，然后乘以剩余部分。 即使对于非大整数，该算法也能很好地工作。
在自然数序列中，第k个数可以整除第(k的倍数)个数。通过利用这个事实，首先取消分母，然后乘以剩余部分。这确保了如果答案没有溢出，在计算过程中也不会溢出。 入山算法 
public static BigInteger Combination(int n, int r)
{
    if (n < 0 || r < 0 || r > n) throw new ArgumentException("Invalid parameter");

    if (n - r < r) r = n - r;
    if (r == 0) return 1;
    if (r == 1) return n;

    int[] numerator = new int[r];
    int[] denominator = new int[r];

    for (int k = 0; k < r; k++)
    {
        numerator[k] = n - r + k + 1;
        denominator[k] = k + 1;
    }

    for (int p = 2; p <= r; p++)
    {
        int pivot = denominator[p - 1];
        if (pivot > 1)
        {
            int offset = (n - r) % p;
            for (int k = p - 1; k < r; k += p)
            {
                numerator[k - offset] /= pivot;
                denominator[k] /= pivot;
            }
        }
    }

    BigInteger result = BigInteger.One;
    for (int k = 0; k < r; k++)
    {
        if (numerator[k] > 1) result *= numerator[k];
    }
    return result;
}
目前你的回答不够清晰。请编辑并添加更多细节，以帮助其他人理解它如何回答所提出的问题。你可以在帮助中心找到有关如何撰写好答案的更多信息。 - Community
加入代码：你现在清楚了吗？ - iriyamanorio
已添加代码：你现在清楚了吗？ - undefined
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接