Расчет суммы комбинаций

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;
}

Этот алгоритм также содержится в книге Кнута «Искусство компьютерного программирования», 3-е издание, том 2: Получисленные алгоритмы ».

ОБНОВЛЕНИЕ: Существует небольшая вероятность того, что алгоритм переполнится в строке:

r *= n--;

для очень большой n. Наивная верхняя граница равна sqrt(std::numeric_limits<long long>::max()), что означает n меньше, чем грубые 4 000 000 000.

Answer 2

От Ответ Андреаса :

Вот древний алгоритм, который точен и не переполняется, если результат не слишком велик для 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--;

для очень большой 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, а затем вычислить г просто:

// 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), вы получите исключение вместо неправильного ответа.

Answer 3

Следующая процедура вычислит 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;
}

Answer 4

Помните, что

n! / ( n - r )! = n * ( n - 1) * .. * (n - r + 1 )

так что оно меньше n !. Таким образом, решение состоит в том, чтобы оценить n * (n - 1) * ... * (n - r + 1) вместо первого вычисления n! а затем разделить его.

Конечно, все зависит от относительной величины n и r - если r относительно велико по сравнению с n, то оно все равно не подходит.

Answer 5

Ну, я должен ответить на свой вопрос. Я читал о треугольнике Паскаля и случайно заметил, что мы можем вычислить количество комбинаций с ним:

#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;
}

Answer 6

Получение простой факторизации биномиального коэффициента, вероятно, является наиболее эффективным способом его вычисления, особенно если умножение стоит дорого. Это, безусловно, относится к связанной с этим проблеме вычисления факториала (см., Например, Нажмите здесь ).

Вот простой алгоритм, основанный на сите Эратосфена, который вычисляет простую факторизацию. Идея в основном состоит в том, чтобы пройти через простые числа, как вы находите их, используя сито, а затем также рассчитать, сколько их кратных попадает в диапазоны [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

Answer 7

Используя грязный трюк с длинным двойником, можно получить ту же точность, что и Говард Хиннант (и, возможно, больше):

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.

Answer 8

Один из кратчайших путей:

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);
}

Answer 9

Если вы хотите быть на 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).

Answer 10

Сначала упростите формулу. Ты не хочешь делать длинное деление.