Остаток от деления больших биномиальных коэффициентов

[malphunction]

Решаю олимпиадную задачку по программированию, для её решения нужно быстро считать модули биномиальных коэффициентов.

Вот исходные данные:
* Нужно подсчитать $C^n_k (mod\ p)$, где
* $0 < k < n < 10^5$
* $0 < p < 10^5$
* Нужно подсчитать 500 коэффициентов за 3 секунды (500 - это количество входных данных), а лучше -- быстрее, т.к. выполняются другие операции
* Для справки: за секунду выполняется $10^9$ сложений, $10^8$ умножений и $10^7$ делений.

Обращаю ваше внимание, что не нужно считать весь коэффициент, нужно посчитать только остаток от деления. Это как-то должно ускорить вычисления, но как??

Пока придумал такой алгоритм: рассмотрим способ вычисления $C^n_k = {{n \times \ldots \times (n - k + 1) } \over {k!}}$. Ясно, что в числителе и знаменателе одинаковое количество чисел; этим и воспользуемся для того, чтобы избавиться от знаменателя. Будем называть числа в числителе "верхними числами", а в знаменателе -- "нижними". Так вот, ясно, что для каждого нижнего числа есть одно или несколько кратных ему верхних.

Если совпадение однозначно, то удается очень быстро избавиться от знаменателя; останется только числитель. Затем нужно быстро пробежаться по нему, чтобы найти модуль этого произведения. Тогда алгоритму требуется два прохода по массиву верхних чисел и $2*10^5$ делений, и $10^5$ умножений, и он укладывается в те самые 3 секунды.

Однако проблема в том, что для некоторых коэффициентов не получается построить однозначное соответствие: во-первых, среди верхних чисел могут быть простые и им не соответствуют никакие нижние; во-вторых, нижние числа могут либо соответствовать нескольким верхним. Поэтому возникает задача перебора (для подбора верхних чисел, кратных нижним), которая решается достаточно долго.

Вот и получается, что по такому алгоритму я могу считать биномиальные коэффициенты в общем случае за 1 секунду, что в 166 раз хуже, чем надо.

Похоже, нужно как-то опираться на тот факт, что вычисления идут по модулю, но как -- не знаю.

Подскажите, пожалуйста!

На всякий случай, исходная задача -- подсчитать количество перестановок $1..N$ с неподвижными точками по модулю $p$, где $0 < N < 10^9$, она сводится к вычислению биномиального коэффициента и субфакториала.

[Хорхе]

Сперва посчитать, какая степень $p$ в числителе, а какая в знаменателе. Если в числителе больше, сразу вывод.

Если одинакова - использовать теорему Вильсона (легко гуглится). Таким образом для каждого биномиального коэффициента останется меньше $2p$ умножений (деления тоже можно заменить умножениями по теореме Вильсона).

Если надо будет, напишу подробней.

Впрочем, я не понял, как эти биномиальные коэффициенты по модулю связаны с перестановками, ну и ладно.

[malphunction]

Хорхе в сообщении #423442 писал(а):

Впрочем, я не понял, как эти биномиальные коэффициенты по модулю связаны с перестановками, ну и ладно.

Сама задачка -- тут: https://www.spoj.pl/problems/SEQN/
Пусть $p_n(k)$ -- количество перестановок размера $n$ с $k$ неподвижными точками. Я смог вывести, что $p_n(k) = C_n^k \times !(n - k)$, где $!n$ -- субфакториал от $n$. Для субфакториала я смог найти эффективное вычисление по модулю, а $C_n^k$ от больших чисел приводится к меньшим по модулю. Засада в том, что даже с маленькими числами они долго считаются.

Про сравнение модулей в числителе/знаменателе и теорему Вильсона -- спасибо, сейчас поищу...

-- Ср мар 16, 2011 16:35:50 --

Хорхе в сообщении #423442 писал(а):

Сперва посчитать, какая степень $p$ в числителе, а какая в знаменателе. Если в числителе больше, сразу вывод.

Если одинакова - использовать теорему Вильсона (легко гуглится). Таким образом для каждого биномиального коэффициента останется меньше $2p$ умножений (деления тоже можно заменить умножениями по теореме Вильсона).

Хм, к сожалению, так и не понял, как это использовать. Поясните, пожалуйста.

Чтобы посчитать степень в числителе/знаменателе, мне нужно будет факторизовать числа оттуда. По-моему, это займет больше времени, чем 3 секунды. Или существуют быстрые способы факторизации?

Нашел только теорему Вильсона про простоту числа, о том, что $n$ простое, если $(n - 1)! + 1 \equiv 0\pmod{n}$. Придумал только один способ применения: разложим $p$ на простые $p_i$, после чего можно будет заменить ${p_{max}}!$ на 1 (по модулю. Сюда же как-то ещё добавить Китайскую теорему об остатках, да?), и это несколько ускорит умножение.

[malphunction]

Понял, как можно сравнительно быстро проверить вхождения $p$ в числитель/знаменатель. Для этого нужно факторизовать одно только $p$ (на простые множители $pp_i$ со степенями $ps_i$). После чего пробежаться по числам числителя/знаменателя и посчитать степени вхождения простых чисел $pp_i$ в эти числа. Получатся степени $up_i$ для числителя и $down_i$ для знаменателя. По этим массивам получаем кратность $pp_i$ и, соответственно, степерь $p$. Так?

Но тогда получается оценка: $10^5$ чисел в числителе/знаменателе, и около 1-17 простых чисел в разложении $p$. Т.е. нужно будет около $10^6 \ldots 10^7$ делений, а это уже на грани допустимого. Т.е. такой способ не очень подходит :(

[Sonic86]

Насколько я помню, $C_n^k \mod p$ можно найти быстро как Хорхе говорит, только там уже была готовая формула, именная, я ее забыл
Что-то вроде $\binom{n}{k} \equiv \binom{n \mod p}{k \mod p} \binom{\left[ \frac{n}{p} \right]}{\left[ \frac{f}{p} \right]} \pmod p$. И доказывалась она так же: представляем $n=pq+r, k=pu+v$, подставляем, выделяем отдельно множители с $pq+1,...,pq+r;pu+1,...,pu+v$ - это дает один биномиальный коэффициент, а во втором с помощью теоремы Вильсона некратные $p$ множители дают 1, а в кратных сокращается $p$ и получается бином. коэффициент, только меньший...

-- Ср мар 16, 2011 16:56:19 --

Вот оно! формула Люка!
http://ru.wikipedia.org/wiki/%D2%E5%EE% ... B%FE%EA%E0

[venco]

Если $p$ простое, то я бы решал так:
Заранее вычислить обратные элементы всех чисел $[1...{10}^5]$ по модулю $p$. Это можно сделать за порядка ${10}^5$ операций.
А дальше вычислять результат, пользуясь умножением на обратный элемент вместо деления.

[malphunction]

Sonic86 в сообщении #423517 писал(а):

Вот оно! формула Люка!
http://ru.wikipedia.org/wiki/%D2%E5%EE% ... B%FE%EA%E0

Эх, я уже использовал её, чтобы привести $C_n^k$ (с $0 < k < n < 10^9$) к произведениям $C_{n'}^{k'}$ (с $0 < k' < n' < p < 10^5$).

Остается только возможность ещё раз использовать эту формулу для множителей $p$, но если $p$ -- простое и большое, то эта формула не помогает.

[venco]

malphunction в сообщении #423451 писал(а):

Но тогда получается оценка: $10^5$ чисел в числителе/знаменателе, и около 1-17 простых чисел в разложении $p$. Т.е. нужно будет около $10^6 \ldots 10^7$ делений, а это уже на грани допустимого.

Зачем так много делений?
Количество вхождений простого делителя в факториал вычисляется по простой формуле.

[Sonic86]

venco писал(а):

Количество вхождений простого делителя в факториал вычисляется по простой формуле.

Угу, вот по такой:

$$\text{ord}_p(n!) = \sum\limits_{k=1}^N \left[ \frac{n}{p^k}\right]$$

$N = [\log _p n]$
Сразу же можно из нее вывести и формулу для $\text{ord}_p(C_n^k)$, если нужно.

[malphunction]

Всем спасибо за ответы! В результате у меня получается следующая схема:

• Разложить $p$ на простые множители $p_i$.
• Для каждого $p_i$ вычислить $c_i = C_n^k \pmod{p_i}$, вычисление такого бином.коэфф. сокращается по теореме Люка к вычислению $C_{n_{j}}^{k_{j}}$, где $n_j < p_i$.
• Всё равно, возможно, что $p_i$ будет большим (или $p$ -- простым), тогда для вычисления модуля б.к. обращаем знаменатель (используя расширенный алгоритм Евклида, например).
• Далее, из $c_i$ сооружаем целиком $C_n^k$, можно для этого использовать метод независимых остатков из "Конкретной математики" (стр.152). Правда, тут проблема в том, что при разложении $p$ могут получаться неединичные степени...

Сложность этого алгоритма не считал, но, кажется, она будет меньше, чем $10^9$ умножений.

Как вы считаете, всё ли правильно?

[Sonic86]

malphunction писал(а):

Всё равно, возможно, что $p_i$ будет большим (или $p$ -- простым), тогда для вычисления модуля б.к. обращаем знаменатель (используя расширенный алгоритм Евклида, например)

Я вот не знаю, если честно, но, по-моему, очевидного преимущества алгоритма Евклида по сравнению с простым делением нету. Обращение числа $j$ по модулю $p$ потребует $O(\ln p)$ операций, для $1 \leq j \leq k$ надо $O(k \ln p)$ операций, а так всего $k$ делений

malphunction писал(а):

Правда, тут проблема в том, что при разложении $p$ могут получаться неединичные степени...

Да - видимо придется сокращать отдельно, а при простом делении - нет.
Но китайская теорема об остатках позволяет сооружать остатки от деления для произведения любых неединичных степеней простых.

Еще мелочь: $C_n^k = C_n^{n-k}$ - можно добиться $k \leq \frac{n}{2}$

-- Вс мар 20, 2011 19:03:35 --

Еще. Пусть мы считаем $C^k_n = \prod\limits_{p_j \leq n} p_j^{a_j}$.
Из формулы для $a_j$ получается, что надо $\frac{\ln n}{\ln p_j}$ операций выполнить, да потом перемножить $O(\pi (n)) = O(\frac{n}{\ln n})$ чисел. Это всего $\sum\limits_{p \leq n}\frac{\ln n}{\ln p_j} + O(\frac{n}{\ln n}) = O(\ln n \frac{n}{\ln n}) + O(\frac{n}{\ln n}) = O(n)$ операций, так что этот способ, видимо медленнее, чем просто $k$ умножений и $k$ делений, хотя точно не уверен...

[malphunction]

Sonic86 в сообщении #425099 писал(а):

Я вот не знаю, если честно, но, по-моему, очевидного преимущества алгоритма Евклида по сравнению с простым делением нету.

Нужно делать вычисления в целочисленной арифметике, оставаясь в пределах машинных слов ($n \in [0..(2^{32}-1)]$ в моём случае). Я пробовал, чисто для проверки, использовать всякие библиотеки длинной арифметики, типа NTL; там можно посчитать б.к. по обычной формуле, разделив громадный числитель на знаменатель, но получается слишком медленно.

Поэтому и хочется как-то исхитриться и остаться в машинной разрядной сетке. Если заменить деление умножением, то мы избавимся от потенциальных вещественных чисел (дробей), и при умножении у нас всегда будут целые числа. Рост произведения можно контролировать, деля его на $p$. То есть в принципе, б.к. можно посчитать где-то за $O(n \log n)$ умножений, но вот коэффициент при этом $O$ будет большой :(

Ну и в самом начале топика я предлагал такой способ: пройтись по числителю/знаменателю и посокращать, сколько можно. Проход требует линейного времени и очень быстр. Но всё равно, часть этого знаменателя остается, где-то $1/50$ чисел. Можно обращать именно эти числа, а не весь знаменатель.

За остальные подсказки -- спасибо, сегодня вечером ещё поразбираюсь...

-- Пн мар 21, 2011 10:25:35 --

venco в сообщении #423580 писал(а):

Заранее вычислить обратные элементы всех чисел $[1...{10}^5]$ по модулю $p$. Это можно сделать за порядка ${10}^5$ операций.

А как можно сделать за ${10}^5$ операций? Если использовать алгоритм бинарного возведения в степень (http://e-maxx.ru/algo/reverse_element, вроде бы он быстрее сходится, чем алг.Евклида), то на каждый элемент потребуется где-то $\log 32$ умножений и делений, т.е. около 10 операций. И самих элементов $10^5$, получается всего $10^6$ операций.

Как же сделать быстрее?

[Sonic86]

malphunction писал(а):

там можно посчитать б.к. по обычной формуле, разделив громадный числитель на знаменатель, но получается слишком медленно.

Ааа, понял! Я туплю, я просто забыл, что $C_n^k$ очень длинный получается...
А вот это вот:

malphunction писал(а):

Ну и в самом начале топика я предлагал такой способ: пройтись по числителю/знаменателю и посокращать, сколько можно. Проход требует линейного времени и очень быстр. Но всё равно, часть этого знаменателя остается, где-то $1/50$ чисел. Можно обращать именно эти числа, а не весь знаменатель.

и аналогичное сообщение в первом посте я не понял. Как это часть знаменателя остается? Она никогда не остается, биномиальный коэффициент-то целый! Это только если Вы на каком-то шаге вычисления не берете остаток от деления - тогда сразу возникает облом.
Но если сокращать и даже игнорировать то, что тогда придется использовать длинные целые, то будет то, что я вот тут написал:

Sonic86 писал(а):

Еще. Пусть мы считаем $C^k_n = \prod\limits_{p_j \leq n} p_j^{a_j}$.... $O(n)$ операций

+ хранить множество простых за $O(\frac{n}{\ln n})$ памяти.

Кстати, я еще вспомнил треугольник Паскаля. $C_n^k = C_{n-1}^k + C_{n-1}^{k-1}$. Никакого деления и умножения, только взятие остатка. Я не знаю, как машинно считается остаток, но если и его считать сложно (через деление), можно его заменить if-ом с одним вычитанием. Только надо с ним похимичить еще наверное...

(Оффтоп)

вообще тема интересная, меня вот интересовало как считать $C_{2n}^n \mod n$, но тогда компа у меня еще не было...

-- Пн мар 21, 2011 13:08:08 --

Наверное я с треугольником Паскаля загнул. Это надо завести массив из $k$ элементов и выполнить в нем $O(kn)=O(n^2)$ операций сложения.
А если попытаться использовать тождество $C_n^k = \sum\limits_{j=0}^k C_k^jC^{k-j}_{n-k}$, то тут будет и умножение, и $O(n)$ сложений, и рекурсивное вычисление...

[malphunction]

Хорхе в сообщении #423442 писал(а):

Если одинакова - использовать теорему Вильсона (легко гуглится). Таким образом для каждого биномиального коэффициента останется меньше $2p$ умножений (деления тоже можно заменить умножениями по теореме Вильсона).

Кто-нибудь, можете объяснить, как применить теорему Вильсона? Особенно для замены деления на умножения??
Хорхе, к сожалению, не отвечает :(

[Sonic86]

malphunction писал(а):

Кто-нибудь, можете объяснить, как применить теорему Вильсона? Особенно для замены деления на умножения??
Хорхе, к сожалению, не отвечает :(

Я могу. Сначала используете формулу Люка и тогда $0 \leq k \leq n <p$, а потом поскольку $(p-1)! \equiv -1 \pmod p$ будет $C_n^k = \frac{n(n-1)...(n-k+1)}{k!} \equiv - \frac{n(n-1)...(n-k+1)(p-1)!}{k!} = $
$= - n(n-1)...(n-k+1) \cdot (k+1)(k+2)...(p-1) \pmod p$

Для модулей $p^a,a>1$ теорема вроде уже не работает