logo标签列表

在C++中计算组合数（N选R）的数量

c++combinatoricsbinomial-coefficients
23

我正在尝试用C++编写一个查找NCR的程序。但是我的结果有问题，不正确。您能帮助我找出程序中的错误吗？

#include <iostream>
using namespace std;
int fact(int n){
    if(n==0) return 1;
    if (n>0) return n*fact(n-1);
};

int NCR(int n,int r){
    if(n==r) return 1;
    if (r==0&&n!=0) return 1;
    else return (n*fact(n-1))/fact(n-1)*fact(n-r);
};

int main(){
    int n;  //cout<<"Enter A Digit for n";
    cin>>n;
    int r;
         //cout<<"Enter A Digit for r";
    cin>>r;
    int result=NCR(n,r);
    cout<<result;
    return 0;
}
3
你的公式有误，在分子和分母中都有 fact(n-1) （它们相互抵消了）。 - Ben Voigt
我投票关闭此问题，因为它似乎根本是一个数学问题，而不是一个编程问题。 - Karl Knechtel
8个回答
40

你的公式完全错误，正确的应该是fact(n)/fact(r)/fact(n-r)，但这种计算方式效率非常低下。

请参见快速计算多分类组合数，尤其是我的评论。(哦，还有，请重新打开那个问题，这样我才能正确回答)

单一拆分情况实际上非常容易处理：

unsigned nChoosek( unsigned n, unsigned k )
{
    if (k > n) return 0;
    if (k * 2 > n) k = n-k;
    if (k == 0) return 1;

    int result = n;
    for( int i = 2; i <= k; ++i ) {
        result *= (n-i+1);
        result /= i;
    }
    return result;
}

演示：http://ideone.com/aDJXNO

如果结果不适合，您可以计算对数的总和，并以double形式获得近似组合数。或者使用任意精度整数库。

我将我的解决方案放在这篇与之密切相关的问题中，因为ideone.com最近丢失了代码片段，而那个问题仍然关闭以防止新答案的提交。

#include <utility>
#include <vector>

std::vector< std::pair<int, int> > factor_table;
void fill_sieve( int n )
{
    factor_table.resize(n+1);
    for( int i = 1; i <= n; ++i )
        factor_table[i] = std::pair<int, int>(i, 1);
    for( int j = 2, j2 = 4; j2 <= n; (j2 += j), (j2 += ++j) ) {
        if (factor_table[j].second == 1) {
            int i = j;
            int ij = j2;
            while (ij <= n) {
                factor_table[ij] = std::pair<int, int>(j, i);
                ++i;
                ij += j;
            }
        }
    }
}

std::vector<unsigned> powers;

template<int dir>
void factor( int num )
{
    while (num != 1) {
        powers[factor_table[num].first] += dir;
        num = factor_table[num].second;
    }
}

template<unsigned N>
void calc_combinations(unsigned (&bin_sizes)[N])
{
    using std::swap;

    powers.resize(0);
    if (N < 2) return;

    unsigned& largest = bin_sizes[0];
    size_t sum = largest;
    for( int bin = 1; bin < N; ++bin ) {
        unsigned& this_bin = bin_sizes[bin];
        sum += this_bin;
        if (this_bin > largest) swap(this_bin, largest);
    }
    fill_sieve(sum);

    powers.resize(sum+1);
    for( unsigned i = largest + 1; i <= sum; ++i ) factor<+1>(i);
    for( unsigned bin = 1; bin < N; ++bin )
        for( unsigned j = 2; j <= bin_sizes[bin]; ++j ) factor<-1>(j);
}

#include <iostream>
#include <cmath>
int main(void)
{
    unsigned bin_sizes[] = { 8, 1, 18, 19, 10, 10, 7, 18, 7, 2, 16, 8, 5, 8, 2, 3, 19, 19, 12, 1, 5, 7, 16, 0, 1, 3, 13, 15, 13, 9, 11, 6, 15, 4, 14, 4, 7, 13, 16, 2, 19, 16, 10, 9, 9, 6, 10, 10, 16, 16 };
    calc_combinations(bin_sizes);
    char* sep = "";
    for( unsigned i = 0; i < powers.size(); ++i ) {
        if (powers[i]) {
            std::cout << sep << i;
            sep = " * ";
            if (powers[i] > 1)
                std::cout << "**" << powers[i];
        }
    }
    std::cout << "\n\n";
}
3
看起来 n = k 的情况将会被处理，方法是：如果 (k * 2 > n)，则 k = n-k; (此时 k 变为0) 并且如果 (k == 0)，则返回1。 - dodecaplex
@IAmRoot永远不会过度复杂化事情。 - taurus05
警告：nChoosek（）是错误的！如果我们在此行之前添加行assert（0 ==（result％i））;，并计算例如nChoosek（64,7），则result / = i;必须可被整除，您将看到此断言失败：53546048 7（即53546048 / 7 = 7649435.428571428），不是整数。 - user873275
Steven的以下NCR()是正确的，它给出了正确的答案：NCR(64, 7) = 621216192。如果在res /= k;之前添加断言assert(0 == (res % k));，则该断言始终成立。 - user873275
@user873275：除非您的输入需要它，否则这是浪费。如果您的输入更大，则使用任意精度算术库可能会更好。但是，如果int或unsigned类型的变量计算不正确，则应扩展变量result的类型和函数的返回类型。 - Ben Voigt
显示剩余7条评论
18

N选R的定义是计算两个乘积，然后用一个除以另一个：

(N * N-1 * N-2 * ... * N-R+1) / (1 * 2 * 3 * ... * R)

然而，乘法可能会变得非常快并且超出现有数据类型的范围。 实现技巧是重新排列乘法和除法，如下所示：

(N)/1 * (N-1)/2 * (N-2)/3 * ... * (N-R+1)/R

每一步都保证结果是可被整除的（对于n个连续数，其中之一必须被n整除，因此这些数字的乘积也是如此）。

例如，对于N选3，至少N、N-1、N-2中的一个将是3的倍数，而对于N选4，至少N、N-1、N-2、N-3中的一个将是4的倍数。

C++代码如下所示。


```cpp long long choose(int n, int k) { long long res = 1; for (int i = 1; i <= k; ++i) { res *= n - i + 1; res /= i; } return res; } ```
int NCR(int n, int r)
{
    if (r == 0) return 1;

    /*
     Extra computation saving for large R,
     using property:
     N choose R = N choose (N-R)
    */
    if (r > n / 2) return NCR(n, n - r); 

    long res = 1; 

    for (int k = 1; k <= r; ++k)
    {
        res *= n - k + 1;
        res /= k;
    }

    return res;
}
res 是否应该是 long 类型？ - ar2015
1
仅适用于在除法发生之前，res *= n - k + 1 步骤溢出的情况。 - Steven
注意：这应该是正确的答案，而不是最高投票的Ben的答案“nChoosek()”。 “NCR（64,7）= 621216192”。如果在“res / = k;”之前添加它，则断言“assert（0 ==（res％k））;”始终保持。 - user873275
@user873275：它执行的操作与我的完全相同。 - Ben Voigt
@user873275：在第一次循环迭代（k == 1）后，Steven得到了什么？哦，是的，它乘以n并除以1，达到了我相同的起点。当然，这个属性成立，Steven花费了他答案的前2/3来解释为什么。 - Ben Voigt
@Ben Voigt：Steven的解决方案使用了res类型long，而您的解决方案使用了result类型int，这会发生溢出现象。 - user873275
3
一种实现n选k的好方法是将其基于“递增乘积”函数，而不是基于阶乘，该函数与阶乘密切相关。
递增乘积(m,n)将m * (m + 1) * (m + 2) * ... * n相乘，在处理各种特殊情况时有规则，例如n> = m或n <= 1：
请参见此处，其中包含一个在C语言中编写的解释型编程语言中作为内在函数的nCk和nPk实现。
static val rising_product(val m, val n)
{
  val acc;

  if (lt(n, one))
    return one;

  if (ge(m, n))
    return one;

  if (lt(m, one))
    m = one;

  acc = m;

  m = plus(m, one);

  while (le(m, n)) {
    acc = mul(acc, m);
    m = plus(m, one);
  }

  return acc;
}

val n_choose_k(val n, val k)
{
  val top = rising_product(plus(minus(n, k), one), n);
  val bottom = rising_product(one, k);
  return trunc(top, bottom);
}

val n_perm_k(val n, val k)
{
  return rising_product(plus(minus(n, k), one), n);
}

这段代码没有使用像+和<这样的运算符，因为它是类型通用的（类型val表示任何类型的值，例如各种类型的数字，包括"bignum"整数），而且它是用C语言编写的（没有重载），因为它是Lisp-like语言的基础，该语言没有中缀语法。尽管如此，这个n-choose-k实现具有简单的结构，容易理解。
说明：le：小于等于；ge：大于等于；trunc：截断除法；plus：加法，mul：乘法，one：一个类型为val的常量，表示数字1。
1

这条线

else return (n*fact(n-1))/fact(n-1)*fact(n-r);

应该是

else return (n*fact(n-1))/(fact(r)*fact(n-r));

或者甚至更好

else return fact(n)/(fact(r)*fact(n-r));
你和原帖作者犯了一些相同的错误... fact(n-r) 需要是一个除数，而不是乘法因子。 - Ben Voigt
1
这是为了在竞赛编程中解决nCr时不超时的参考，我发布这篇文章是因为它对你有帮助，因为你已经得到了问题的答案。获取二项式系数的质因数分解可能是计算它最有效的方法，特别是如果乘法很昂贵。这在计算阶乘的相关问题上也是正确的（例如，请参见Click here）。以下是一个基于埃拉托色尼筛法的简单算法，用于计算质因数分解。基本思想是通过使用筛法找到素数并计算其倍数在范围[1，k]和[n-k+1，n]中的数量。Sieve本质上是一个O(nloglogn)的算法，但没有进行乘法。一旦发现质因数分解，实际所需的乘法次数最坏情况下为O（nloglogn / logn），可能还有更快的方法。
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
请注意，我在三年前的答案中提供的代码不仅可以在不使用动态列表增长的情况下完成此操作，而且还使用因式分解来解决组合数问题。而且它是用问题所要求的语言编写的（我猜这里是Python？）。 - Ben Voigt
1

使用double代替int。

更新：

您的公式也是错误的。您应该使用fact(n)/fact(r)/fact(n-r)。

是的！我没有在意公式，因为使用“int”也会在那里产生错误的结果。 - Pulkit Goyal
3
即使使用“double”也不是一个很好的方法。尝试计算“6000选择3”。它很容易适应32位“int”，但使用该公式和双倍将失败。 - Ben Voigt
哎呀，我的意思是 600 中选 3 很容易。但是 6000 中选 3 实际上不行。 - Ben Voigt
0
// CPP program To calculate The Value Of nCr
#include <bits/stdc++.h>
using namespace std;

int fact(int n);

int nCr(int n, int r)
{
    return fact(n) / (fact(r) * fact(n - r));
}

// Returns factorial of n
int fact(int n)
{
    int res = 1;
    for (int i = 2; i <= n; i++)
        res = res * i;
    return res;
}

// Driver code
int main()
{
    int n = 5, r = 3;
    cout << nCr(n, r);
    return 0;
}
0

递归函数在这里使用不正确。fact() 函数应该改成这样：

int fact(int n){
if(n==0||n==1) //factorial of both 0 and 1 is 1. Base case.
{
    return 1;
}else

    return (n*fact(n-1));//recursive call.

};

递归调用应该在else部分进行。

NCR()函数应该改为：

int NCR(int n,int r){
    if(n==r) {
        return 1;
    } else if (r==0&&n!=0) {
        return 1;
    } else if(r==1)
    {
        return n;
    }
    else
    {
        return fact(n)/(fact(r)*fact(n-r));
    }
};
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接