计算组合的数量

Question

干杯，

我知道你可以用下面的公式得到组合的数量（不重复，顺序并不重要）：

// Choose r from n

n! / r!(n - r)!

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

n = 52

n! = 8,0658175170943878571660636856404e+67

数字变得太大了对于unsigned __int64（或unsigned long long）。是否有一些解决方法可以在没有任何第三方“bigint”库的情况下实现该公式？

Answer 1

这是一个古老的算法，它是精确的，并且不会溢出，除非结果对于 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::max())，这意味着 n 大约小于 4,000,000,000。

Answer 2

来自Andreas的回答：

这是一种古老的算法，它非常精确并且不会溢出，除非结果对于 long long 来说太大

无符号长长
选择（无符号长长n，无符号长长k）{
    如果（k>n）{
        返回0；
    }
    无符号长长 r = 1;
    for (无符号长长 d = 1; d <= k; ++d) {
        r *= n--;
        r/=d；
    }
    返回 r；
}

我认为这个算法也在 Knuth 的《计算机编程艺术，第 3 版，第 2 卷：半数值算法》中。

更新：算法在线上溢出的可能性很小：

<前><代码>r *= n--;

对于非常大的n。简单的上限是 sqrt(std::numeric_limits::max()) ，这意味着 n 大约小于 4,000,000,000。

考虑 n == 67 和 k == 33。上述算法会因 64 位 unsigned long long 溢出。然而，正确答案可以用 64 位表示：14,226,520,737,620,288,370。上述算法对其溢出保持沉默，选择（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，因此我们不需要计算它。因此我们可以将 n 除以 d_temp：

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

现在您可以计算 select(67, 33) 而不会溢出。如果您尝试选择(68, 33)，您将得到一个异常而不是错误的答案。

Answer 3

以下例程将使用递归定义和记忆来计算 n-choose-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;
}

Answer 4

请记住

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

所以它比 n! 小得多。所以解决方案是评估 n* ( n - 1 ) * ... * ( n - r + 1) 而不是先计算 n！然后划分它。

当然，这完全取决于 n 和 r 的相对大小 - 如果 r 与 n 相比相对较大，那么它仍然不适合。

Answer 5

好吧，我必须回答我自己的问题。我在阅读有关帕斯卡三角形的内容时，偶然注意到我们可以用它计算组合的数量：

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

Answer 6

稍微改进一下 Howard Hinnant 的答案（在这个问题中）：
每个循环调用 gcd() 似乎有点慢。
我们可以将 gcd() 调用聚合到最后一个，同时充分利用 Knuth 的书“计算机编程的艺术，第 3 版，第 2 卷：半数值算法”中的标准算法：

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() 3次。

如果您愿意，请将 uint64_t 替换为 unsigned long long。

Answer 7

获得二项式系数的素因数分解可能是计算它的最有效的方法，特别是在乘法成本昂贵的情况下。这对于计算阶乘的相关问题来说确实如此（例如，请参见单击此处 ）。

这是一个基于埃拉托斯特尼筛法的简单算法，用于计算素因数分解。这个想法基本上是在使用筛子找到质数时遍历它们，然后计算它们的倍数中有多少落在 [1, k] 和 [n-k+1,n] 范围内。筛子本质上是一个 O(n \log \log n) 算法，但没有进行乘法运算。一旦找到质因数分解，所需的实际乘法次数最多是 O\left(\frac{n \log \log n}{\log n}\right) 并且可能有比这更快的方法。

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

Answer 8

使用带有 long double 的肮脏技巧，可以获得与 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 进行近似。

Answer 9

最短的方法之一：

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

Answer 10

类似于埃拉托斯特尼筛法的方法。
虽然埃拉托斯特尼的筛子是多重湮灭，
这一次是多重半杀。
由于 n!/((nr)!r!) 始终是整数，因此首先取消分母，然后乘以其余部分。
即使对于非大整数，该算法也能很好地工作。

在自然数序列中，第k个数可以整除第（k的倍数）个数。这可以通过 k=2,3,4,... 连续完成，利用这个事实，首先取消分母，然后乘以余数。这确保了如果答案没有溢出，则在计算过程中也不会溢出。

入山算法

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

Answer 11

如果你想百分百确定只要最终结果在数值限制内就不会发生溢出，你可以逐行求帕斯卡三角求和：

for (int i=0; i<n; i++) {
    for (int j=0; j<=i; j++) {
        if (j == 0) current_row[j] = 1;
        else current_row[j] = prev_row[j] + prev_row[j-1];
    }
    prev_row = current_row; // assume they are vectors
}
// result is now in current_row[r-1]

但是，这种算法比乘法慢得多。因此，也许您可​​以使用乘法来生成所有您知道“安全”的情况，然后从那里使用加法。 （..或者你可以只使用 BigInt 库）。