乾杯、
次の式で組み合わせの量を取得できることはわかっています（繰り返しなしで、順序は重要ではありません）：
// n から r を選ぶ ん！/ r!(n - r)!
ただし、これを C++ で実装する方法がわかりません。
n = 52 ん！= 8,0658175170943878571660636856404e+67
unsigned __int64(または)に対しても数が大きくなりすぎunsigned long longます。サードパーティの「bigint」ライブラリなしで式を実装するための回避策はありますか?
これは、結果が大きすぎない限り、正確でオーバーフローしない古代のアルゴリズムです。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 の「The Art of Computer Programming, 3rd Edition, Volume 2: Seminumerical Algorithms」にもあると思います。
更新:アルゴリズムが回線上でオーバーフローする可能性がわずかにあります。
r *= n--;
非常に大きなn の場合。ナイーブな上限は、およそ 4,000,000,000 未満でsqrt(std::numeric_limits<long long>::max())あることを意味します。n
アンドレアスの答えから：
これは、結果が大きすぎない限り、正確でオーバーフローしない古代のアルゴリズムです。
long longunsigned 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 の「The Art of Computer Programming, 3rd Edition, Volume 2: Seminumerical Algorithms」にもあると思います。
更新:アルゴリズムが回線上でオーバーフローする可能性がわずかにあります。
r *= n--;非常に大きなn の場合。ナイーブな上限は、およそ 4,000,000,000 未満で
sqrt(std::numeric_limits<long long>::max())あることを意味します。n
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* ( n - 1 ) * ... * ( n - r + 1) を評価することです! そしてそれを分割します。
もちろん、それはすべて 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;
}二項係数の素因数分解を取得することは、特に乗算が高価な場合、おそらく最も効率的な計算方法です。これは、階乗の計算に関連する問題にも当てはまります (たとえば、ここをクリックしてください)。
これは、素因数分解を計算するエラトステネスのふるいに基づく単純なアルゴリズムです。基本的には、ふるいを使用して見つけた素数を調べるだけでなく、[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 の回答 (in this question) を少し改善します: ループごとに gcd() を呼び出すのは少し遅いようです。Knuth の著書「The Art of Computer Programming, 3rd Edition, Volume 2: Seminumerical Algorithms」の標準アルゴリズムを最大限に活用しながら、gcd() 呼び出しを最後の呼び出しに集約できます。
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 になることは実際には見られません。組み合わせによる削減はかなり適切です)。
必要に応じて、uint64_t を unsigned long long に置き換えます。
最短の方法の 1 つ:
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% 確実にしたい場合は、パスカルの三角形を行ごとに合計できます。
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 ライブラリを使用することもできます)。