Приветствия,
Я знаю, что вы можете получить количество комбинаций со следующей формулой (без повторения и порядка не важно):
// Choose r from n n! / r!(n - r)!
Однако я не знаю, как реализовать это на С++, поскольку, например, с помощью
n = 52 n! = 8,0658175170943878571660636856404e+67
число становится слишком большим даже для unsigned __int64 (или unsigned long long). Есть ли способ обхода формулы без каких-либо сторонних "bigint" -libraries?
Здесь древний алгоритм, который является точным и не переполняется, если результат не является большим для a 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.
Из ответ Андреаса:
Здесь древний алгоритм, который является точным и не переполняется, если результат не является большим для a
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 "Искусство программирования, 3-е издание, том 2: Семинумерные алгоритмы". Я думаю.
ОБНОВЛЕНИЕ:. Небольшая вероятность того, что алгоритм переполнится в строке:
r *= n--;для очень большого n. Наивная верхняя граница
sqrt(std::numeric_limits<long long>::max()), что означает, чтоnменьше, чем на 4 000 000 000.
Рассмотрим n == 67 и k == 33. Вышеупомянутый алгоритм переполняется 64-битным беззнаковым длинным длинным. И все же правильный ответ представляется в 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, а затем вычислять of 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;
}
Теперь вы можете вычислить select (67, 33) без переполнения. И если вы попытаетесь выбрать (68, 33), вы получите исключение вместо неправильного ответа.
Помните, что
n! / ( n - r )! = n * ( n - 1) * .. * (n - r + 1 )
поэтому он меньше n!. Таким образом, решение состоит в том, чтобы сначала оценить n * (n - 1) *... * (n - r + 1) вместо первого вычисления n! и затем разделив его.
Конечно, все зависит от относительной величины n и r - если r относительно велико по сравнению с n, то это все равно не подойдет.
Следующая процедура будет вычислять n-select-k, используя рекурсивное определение и memoization. Подпрограмма чрезвычайно быстро и точно:
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;
}
Ну, я должен ответить на свой вопрос. Я читал о треугольнике Паскаля и случайно заметил, что мы можем рассчитать количество комбинаций с ним:
#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 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]. Сито по существу является алгоритмом 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
Один из способов SHORTEST:
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).