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

Ответ 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.

Это должно сделать его проще, даже если 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;
}

Ответ 2

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

Обратите внимание, что если 1009 (включая его кратность), появляется больше цифр в числителе, чем знаменатель, тогда mod 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).

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

Наверху мы имеем 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!

Ответ 3

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

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

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

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

(где 0 <= mi, ni < 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

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

Ответ 4

psudo-код для вычисления nCk:

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

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

Цикл идет 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, но этот алгоритм вообще не является тривиальным.

Ответ 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

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