Number of combinations (N choose R) in C++
Here I try to write a program in C++ to find NCR. But I've got a problem in the result. It is not correct. Can you help me find what the mistake is in the program?
#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;
}
Solution 1:
Your formula is totally wrong, it's supposed to be fact(n)/fact(r)/fact(n-r)
, but that is in turn a very inefficient way to compute it.
See Fast computation of multi-category number of combinations and especially my comments on that question. (Oh, and please reopen that question also so I can answer it properly)
The single-split case is actually very easy to handle:
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;
}
Demo: http://ideone.com/aDJXNO
If the result doesn't fit, you can calculate the sum of logarithms and get the number of combinations inexactly as a double. Or use an arbitrary-precision integer library.
I'm putting my solution to the other, closely related question here, because ideone.com has been losing code snippets lately, and the other question is still closed to new answers.
#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";
}
Solution 2:
The definition of N choose R is to compute the two products and divide one with the other,
(N * N-1 * N-2 * ... * N-R+1) / (1 * 2 * 3 * ... * R)
However, the multiplications may become too large really quick and overflow existing data type. The implementation trick is to reorder the multiplication and divisions as,
(N)/1 * (N-1)/2 * (N-2)/3 * ... * (N-R+1)/R
It's guaranteed that at each step the results is divisible (for n continuous numbers, one of them must be divisible by n, so is the product of these numbers).
For example, for N choose 3, at least one of the N, N-1, N-2 will be a multiple of 3, and for N choose 4, at least one of N, N-1, N-2, N-3 will be a multiple of 4.
C++ code given below.
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;
}