干杯，
我知道你可以用下面的公式得到组合的数量（不重复，顺序不重要）：
// 从 n 中选择 r 嗯！ /r!(n - r)!
但是，我不知道如何在 C++ 中实现这一点，例如使用
n = 52 嗯！ = 8,0658175170943878571660636856404e+67
即使对于unsigned __int64（或unsigned long long），这个数字也太大了。是否有一些解决方法可以在没有任何第三方“bigint”-库的情况下实现公式？
【问题讨论】：
标签： c++ algorithm combinatorics
这是一个古老的算法，它是精确的并且不会溢出，除非结果对于 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改进吗，先做除法？
for very large n. A naive upper bound... 开始
这是一个古老的算法，它是精确的并且不会溢出，除非结果对于
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。
考虑 n == 67 和 k == 33。上述算法以 64 位 unsigned long long 溢出。然而，正确答案可以用 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，因此我们不需要计算它。所以我们可以将 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;
}
现在您可以计算choose(67, 33) 而不会溢出。如果你尝试choose(68, 33)，你会得到一个异常而不是错误的答案。
【讨论】：
以下例程将使用递归定义和记忆来计算 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;
}
【讨论】：
记住
n! / ( n - r )! = n * ( n - 1) * .. * (n - r + 1 )
所以它比 n! 小得多。所以解决方案是评估 n* ( n - 1 ) * ... * ( n - r + 1) 而不是首先计算 n！然后分割它。
当然这完全取决于 n 和 r 的相对大小 - 如果 r 与 n 相比相对较大，那么它仍然不适合。
【讨论】：
好吧，我必须回答我自己的问题。我正在阅读有关帕斯卡三角形的内容，偶然注意到我们可以用它计算组合的数量：
#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;
}【讨论】：
uint64_t 是 #include <cstdint> 的一部分，因此我们不再需要在此示例中使用 boost
获取二项式系数的素数分解可能是计算它的最有效方法，尤其是在乘法代价高昂的情况下。计算阶乘的相关问题当然是正确的（例如参见Click here）。
这是一个基于埃拉托色尼筛法的简单算法，用于计算素数分解。这个想法基本上是在使用筛子找到素数时遍历它们，然后还要计算它们有多少倍数落在 [1, k] 和 [n-k+1,n] 范围内。 Sieve 本质上是一个 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
【讨论】：
使用长双打的肮脏技巧，可以获得与 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 进行逼近。
【讨论】：
稍微改进了 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。
【讨论】：
r x b / c
最短的方法之一：
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);
}
【讨论】：
如果您想 100% 确定只要最终结果在数值限制内就不会发生溢出，您可以逐行总结 Pascal 三角：
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 库）。
【讨论】：
n-- 的过程中可能会发生溢出。它不会在这里发生。