计算组合的数量

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<long long>::max())意味着n小于4,000,000,000的粗糙.

Answer 2

来自Andreas的回答:

这是一个古老的算法,它是精确的,并且不会溢出,除非结果是大的 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位无符号长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,因此我们不需要计算它.所以我们可以用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;
}

现在你可以计算选择(67,33)而不会溢出.如果你尝试选择(68,33),你会得到一个例外,而不是一个错误的答案.

Answer 3

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