Вычисление биномиального коэффициента (nCk) для больших n & k
Я только видел этот вопрос и понятия не имею, как его решить. Можете ли вы предоставить мне алгоритмы, коды C ++ или идеи?
Это очень простая проблема. Учитывая значения N и K, нужно сообщить нам значение биномиального коэффициента C (N, K). Вы можете быть уверены, что K < = N и максимальное значение N составляет 1 000 000 000 000 000. Поскольку значение может быть очень большим, вам нужно вычислить результат по модулю 1009. Входные данные
Первая строка входных данных содержит количество тестовых примеров T, не более 1000. Каждая из следующих T строк состоит из двух целых чисел, разделенных пробелами N и K, где 0 <= K <= N и 1 <= N <= 1 000 000 000 000 000. Выходные данные
Для каждого тестового примера выведите на новой строке значение биномиального коэффициента C (N, K) по модулю 1009. Пример
Входной сигнал:
3
3 1
5 2
10 3Выход:
3
10
120
4 ответа
Обратите внимание, что 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 select b, где a <= 1009 и b <= 1009.
Это должно упростить задачу, даже если N близко к 10 ^ 15.
Примечание:
Для N = 10 ^ 15, N выберите N / 2 больше, чем 2 ^ (100000000000000) который путь за беззнаковый длинный длинный.
Кроме того, алгоритм, предложенный Теорема Лукаса равна O (log N), что является
экспоненциальнобыстрее, чем пытались вычислить биномиальный коэффициент напрямую (даже если вы сделали мод 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;
}
Биномиальный коэффициент - это один факториал, деленный на два других, хотя k! член внизу сокращается очевидным образом.
Обратите внимание, что если 1009 (включая его кратные) встречается в числителе больше раз, чем в знаменателе, то ответ по модулю 1009 равен 0. Оно не может появляться в знаменателе больше раз, чем в числителе (поскольку биномиальные коэффициенты являются целыми числами), поэтому единственные случаи, когда вам нужно что-то делать, - это когда оно встречается одинаковое количество раз в обоих. Не забывайте считать, что кратные (1009) ^ 2 равны двум и т. Д.
После этого, я думаю, вы просто убираете мелкие случаи (то есть небольшое количество значений для умножения / деления), хотя я не уверен без нескольких тестов. Положительная сторона 1009 - простое число, поэтому арифметика по модулю 1009 выполняется в поле, что означает, что после исключения кратных 1009 сверху и снизу, вы можете выполнить остальную часть умножения и деления по модулю 1009 в любом порядке.
Там, где остались немаленькие случаи, они все равно будут включать в себя перемножение длинных серий последовательных целых чисел. Это можно упростить, зная 1008! (мод. 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!
Я не думаю, что вы хотите вычислить C (n, k), а затем уменьшить mod 1009. Для самого большого, C (1e15,5e14) потребуется что-то вроде 1e16 бит ~ 1000 терабайт.
Более того, выполнение цикла в ответе snakiles 1e15 раз может занять некоторое время. Вы можете использовать if
n = n0 + n1 * p + n2 * p ^ 2 ... + nd * p ^ d
m = m0 + m1 * p + m2 * p ^ 2 .. . + md * p ^ d
(где 0 <= mi, ni , тогда
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.
псудокод для вычисления nCk:
result = 1
for i=1 to min{K,N-K}:
result *= N-i+1
result /= i
return result
Сложность времени: O (min {K, NK})
Цикл идет от i = 1 до min {K, NK } вместо от i = 1 до K , и это нормально, потому что
C(k,n) = C(k, n-k)
И вы можете вычислить это еще более эффективно, если используете функцию GammaLn.
nCk = exp(GammaLn(n+1)-GammaLn(k+1)-GammaLn(n-k+1))
Функция GammaLn является натуральным логарифмом гамма-функции . Я знаю, что есть эффективный алгоритм для вычисления функции GammaLn, но этот алгоритм совсем нетривиален.