Расчет биномиального коэффициента (nCk) для больших n & k

Question

Я только что видел этот вопрос и понятия не имею, как его решить. Можете ли вы предоставить мне алгоритмы, коды C ++ или идеи?

Это очень простая проблема. Учитывая значения N и K, вам необходимо сообщить нам значение биномиального коэффициента C (N, K). Вы можете быть уверены, что K <= N и максимальное значение N составляет 1 000 000 000 000 000. Поскольку значение может быть очень большим, вам нужно вычислить результат по модулю 1009. Входной </p>

В первой строке входных данных содержится количество тестовых примеров T, не более 1000. Каждая из следующих T строк состоит из двух целых чисел, разделенных пробелами N и K, где 0 <= K <= N и 1 <= N <= 1 000 000 000 000 000 Выход </p>

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

Вход:
3
3 1
5 2
10 3

Выход:
3
10
120

Answer 1

Обратите внимание, что 1009 - простое число.

Теперь вы можете использовать Теорема Лукаса .

Который гласит:

Let p be a prime.
If n  = a1a2...ar when written in base p and
if k  = b1b2...br when written in base p

(pad with zeroes if required)

Then

(n choose k) modulo p = (a1 choose b1) * (a2 choose  b2) * ... * (ar choose br) modulo p.

i.e. remainder of n choose k when divided by p is same as the remainder of
the product (a1 choose b1) * .... * (ar choose br) when divided by p.
Note: if bi > ai then ai choose bi is 0.

Таким образом, ваша задача сводится к тому, чтобы найти произведение по модулю 1009 не более чем log n / log 1009 чисел (количество цифр N в базе 1009) формы a, выберите b, где a <= 1009 и b <= 1009. </p>

Это должно облегчить работу, даже когда N близко к 10 ^ 15.

Примечание:

Для N = 10 ^ 15, N выбрать N / 2 больше 2 ^ (100000000000000), который является способом за неподписанным длинным длинным.

Кроме того, алгоритм, предложенный Теорема Лукаса есть O (log N), которая exponentially быстрее, чем пытаться вычислить биномиальный коэффициент напрямую (даже если вы сделали мод 1009 позаботиться о переполнении).

Вот некоторый код для Binomial, который я написал давно, все, что вам нужно сделать, это изменить его для выполнения операций по модулю 1009 (могут быть ошибки и не обязательно рекомендуемый стиль кодирования):

class Binomial
{
public:
    Binomial(int Max)
    {
        max = Max+1;
        table = new unsigned int * [max]();
        for (int i=0; i < max; i++)
        {
            table[i] = new unsigned int[max]();

            for (int j = 0; j < max; j++)
            {
                table[i][j] = 0;
            }
        }
    }

    ~Binomial()
    {
        for (int i =0; i < max; i++)
        {
            delete table[i];
        }
        delete table;
    }

    unsigned int Choose(unsigned int n, unsigned int k);

private:
    bool Contains(unsigned int n, unsigned int k);

    int max;
    unsigned int **table;
};

unsigned int Binomial::Choose(unsigned int n, unsigned int k)
{
    if (n < k) return 0;
    if (k == 0 || n==1 ) return 1;
    if (n==2 && k==1) return 2;
    if (n==2 && k==2) return 1;
    if (n==k) return 1;


    if (Contains(n,k))
    {
        return table[n][k];
    }
    table[n][k] = Choose(n-1,k) + Choose(n-1,k-1);
    return table[n][k];
}

bool Binomial::Contains(unsigned int n, unsigned int k)
{
    if (table[n][k] == 0) 
    {
        return false;
    }
    return true;
}

Answer 2

Биномиальный коэффициент - это один факториал, деленный на два других, хотя термин k! внизу очевидным образом отменяет.

Обратите внимание, что если 1009 (включая его кратные), появляется больше раз вв числителе, чем в знаменателе, то ответный мод 1009 равен 0. Он не может появляться в знаменателе больше раз, чем в числителе (поскольку биномиальные коэффициенты являются целыми числами), поэтому единственные случаи, когда вы должны что-либо делать, это когдаодинаковое количество раз в обоих.Не забудьте посчитать кратные (1009) ^ 2 как два, и т. Д.

После этого, я думаю, вы просто зачищаете небольшие дела (имеется в виду небольшое количество значений для умножения / деления)Хотя я не уверен без нескольких тестов.На плюсовой стороне 1009 простое число, поэтому в поле происходит арифметическое по модулю 1009, что означает, что после раздачи кратных 1009 сверху и снизу вы можете выполнить остальную часть модуля умножения и деления 1009 в любом порядке.

Там, где остаются немалые дела, они все равно будут включать умножение длинных серий последовательных чисел.Это можно упростить, зная 1008! (mod 1009).Это -1 (1008, если хотите), поскольку 1 ... 1008 - это p-1 ненулевые элементы простого поля над p.Поэтому они состоят из 1, -1, а затем (p-3)/2 пар мультипликативных инверсий.

Так, например, рассмотрим случай C ((1009 ^ 3), 200).

Представьте себечто число 1009 равно (не знаю, если они, потому что я не закодировал формулу, чтобы выяснить), так что это случай, требующий работы.

На вершине у нас есть 201... 1008, который мы должны рассчитать или найти в предварительно вычисленной таблице, затем 1009, затем 1010 ... 2017, 2018, 2019 ... 3026, 3027 и т. Д. Все диапазоны ...1, поэтому нам просто нужно знать, сколько таких диапазонов.

Это оставляет 1009, 2018, 3027, которые, как только мы отменили их со 1009 снизу, будут просто 1, 2, 3, ... 1008, 1010, ..., плюс несколько кратных1009 ^ 2, который мы снова отменим и оставим себя с последовательными целыми числами для умножения.

Мы можем сделать что-то очень похожее с нижней частью, чтобы вычислить мод продукта 1009 из "1 ... 1009 ^ 3 -200 со всеми полномочиями 1009 разделены ".Это оставляет нас с разделением в основной области.IIRC в принципе сложно, но 1009 - это достаточно малое число, чтобы мы могли управлять 1000 из них (верхний предел количества тестовых случаев).

Конечно, при k = 200, существует огромное перекрытие, котороеможет быть отменено более напрямую.Это то, что я имел в виду под маленькими и немалыми случаями: я рассматривал это как немалый случай, когда на самом деле мы могли бы просто уйти от этого "грубым" путем вычисления ((1009^3-199) * ... * 1009^3) / 200!

Answer 3

Не думаю, что вы хотите вычислить C (n, k), а затем уменьшить мод 1009. Самый большой, C (1e15,5e14) потребует что-то вроде 1e16 битов ~ 1000 терабайт

Более того, выполнение цикла в ответе snakiles 1e15 раз может занять некоторое время. Что вы могли бы использовать, если

n = n0 + n1 * p + n2 * p ^ 2 ... + nd * p ^ d

m = m0 + m1 * p + m2 * p ^ 2 ... + md * p ^ d

(где 0 <= mi, ni <p) </p>

, то C (n, m) = C (n0, m0) * C (n1, m1) * ... * C (nd, nd) mod p

см., Например, http://www.cecm.sfu.ca/organics/papers/granville/paper/binomial/html/binomial.html

Одним из способов было бы использование треугольника Паскаля для построения таблицы всех значений C (m, n) для 0 <= m <= n <= 1009. </p>

Answer 4

псевдокод для расчета nCk:

result = 1    
for i=1 to min{K,N-K}:
   result *= N-i+1
   result /= i
return result

Сложность времени: O (min {K, NK})

Цикл идет from i=1 to min{K,N-K} вместо from i=1 to K, и этохорошо, потому что

C(k,n) = C(k, n-k)

И вы можете вычислить эту вещь еще эффективнее, если будете использовать функцию GammaLn.

nCk = exp(GammaLn(n+1)-GammaLn(k+1)-GammaLn(n-k+1))

Функция GammaLn является натуральным логарифмом функции Gamma.Я знаю, что есть эффективный алгоритм для вычисления функции GammaLn, но этот алгоритм совсем не тривиален.

Answer 5

Следующий код показывает, как получить все биномиальные коэффициенты для заданного размера 'n'.Вы можете легко изменить его, чтобы он остановился на данном k, чтобы определить nCk.Он очень эффективен в вычислительном отношении, прост в кодировании и работает для очень больших n и k.

binomial_coefficient = 1
output(binomial_coefficient)
col = 0
n = 5

do while col < n
    binomial_coefficient = binomial_coefficient * (n + 1 - (col + 1)) / (col + 1)
    output(binomial_coefficient)
    col = col + 1
loop

Поэтому вывод биномиальных коэффициентов:

1
1 *  (5 + 1 - (0 + 1)) / (0 + 1) = 5 
5 *  (5 + 1 - (1 + 1)) / (1 + 1) = 15
15 * (5 + 1 - (2 + 1)) / (2 + 1) = 15 
15 * (5 + 1 - (3 + 1)) / (3 + 1) = 5 
5 *  (5 + 1 - (4 + 1)) / (4 + 1) = 1

Я нашелформула когда-то давно в Википедии, но по какой-то причине ее уже нет: (