Расчет суммы комбинаций
Ура,
Я знаю, что вы можете получить количество комбинаций по следующей формуле (без повторения и порядка это не важно):
// Выберите r из n п! / r!(n - r)!
Тем не менее, я не знаю, как реализовать это в C++, так как, например, с
n = 52 п! = 8,0658175170943878571660636856404e+67
число становится слишком большим даже для unsigned __int64 (или же unsigned long long). Есть ли обходной путь для реализации формулы без каких-либо сторонних библиотек "bigint"?
12 ответов
Вот древний алгоритм, который является точным и не переполняется, если результат не является большим для 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;
}
Я думаю, что этот алгоритм также содержится в книге Кнута "Искусство компьютерного программирования", 3-е издание, том 2: Полу численные алгоритмы.
ОБНОВЛЕНИЕ: есть небольшая вероятность того, что алгоритм переполнится на линии:
r *= n--;
для очень большого п. Наивная верхняя граница sqrt(std::numeric_limits<long long>::max()) что означает n менее чем 4 000 000 000
Из ответа Андреаса:
Вот древний алгоритм, который является точным и не переполняется, если результат не является большим для
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; }Я думаю, что этот алгоритм также содержится в книге Кнута "Искусство компьютерного программирования", 3-е издание, том 2: Полу численные алгоритмы.
ОБНОВЛЕНИЕ: есть небольшая вероятность того, что алгоритм переполнится на линии:
r *= 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. И вышеприведенный алгоритм ничего не говорит о его переполнении, select (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;
}
Теперь вы можете вычислить select (67, 33) без переполнения. И если вы попытаетесь выбрать (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 - 1) * ... * ( n - r + 1) вместо первого вычисления n! а затем разделить его.
Конечно, все зависит от относительной величины n и r - если r относительно велико по сравнению с n, то оно все равно не подходит.
Немного улучшает ответ Ховарда Хиннанта (в этом вопросе): вызов gcd() для каждого цикла кажется немного медленным. Мы могли бы объединить вызов gcd() с последним вызовом, максимально используя стандартный алгоритм из книги Кнута «Искусство компьютерного программирования, 3-е издание, том 2: получисловые алгоритмы»:
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, комбинаторное сокращение вполне прилично.), Т.е. вызов select () 3 раза, для случаев без переполнения.
Замените uint64_t на unsigned long long, если хотите.
Ну, я должен ответить на свой вопрос. Я читал о треугольнике Паскаля и случайно заметил, что мы можем рассчитать количество комбинаций с ним:
#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]. Сито по сути является алгоритмом 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
Используя грязный трюк с длинным двойником, можно получить ту же точность, что и Говард Хиннант (и, возможно, больше):
unsigned long long n_choose_k(int n, int k)
{
long double f = n;
for (int i = 1; i<k+1; i++)
f /= i;
for (int i=1; i<k; i++)
f *= n - i;
unsigned long long f_2 = std::round(f);
return f_2;
}
Идея состоит в том, чтобы сначала разделить на k! а затем умножить на n(n-1)...(n-k+1). Аппроксимации через двойник можно избежать, инвертировав порядок цикла for.
Метод, аналогичный Решету Эратосфена. Если решето Эратосфена — это многократное уничтожение, то это — многократное полуубийство. Поскольку n!/((nr)!r!) всегда целое число, сначала сократим знаменатель, а затем умножим остальное. Этот алгоритм хорошо работает даже для небольших целых чисел.
В последовательности натуральных чисел k-е число может делить (кратное k)-е число. Это можно делать непрерывно при k=2,3,4,... Воспользовавшись этим фактом, сначала сократим знаменатель, а затем умножим остаток. Это гарантирует, что если ответ не переполнится, то он не переполнится и в процессе расчета.
public static BigInteger Combination(int n, int r)
{
if (n < 0 || r < 0 || r > n) throw new ArgumentException("Invalid parameter");
if (n - r < r) r = n - r;
if (r == 0) return 1;
if (r == 1) return n;
int[] numerator = new int[r];
int[] denominator = new int[r];
for (int k = 0; k < r; k++)
{
numerator[k] = n - r + k + 1;
denominator[k] = k + 1;
}
for (int p = 2; p <= r; p++)
{
int pivot = denominator[p - 1];
if (pivot > 1)
{
int offset = (n - r) % p;
for (int k = p - 1; k < r; k += p)
{
numerator[k - offset] /= pivot;
denominator[k] /= pivot;
}
}
}
BigInteger result = BigInteger.One;
for (int k = 0; k < r; k++)
{
if (numerator[k] > 1) result *= numerator[k];
}
return result;
}
Сначала упростите формулу. Ты не хочешь делать длинные деления.
Один из кратчайших путей:
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).