Мощное свойство биномиальных коэффициентов, или как сдать задачу на формулу, если последнюю не хочется выводить?
Difference between ru6 and ru7, changed 1,036 character(s)
Привет всем! ↵

Года два назад я решал задачу, в которой для сдачи требовалось вывести формулу, сворачивая суммы, ряды и делая кучу неприятных дел, на которые уходит много сил и энергии на контесте. Я тогда заинтересовался, а существует ли более простой способ вывести эту формулу, не делая почти никаких математических преобразований?↵

[cut]↵

Задолго до этого я наткнулся на [одно свойство биномиальных коэффициентов](https://en.wikipedia.org/wiki/Binomial_coefficient#Binomial_coefficients_as_a_basis_for_the_space_of_polynomials) в Википедии. Я естественно подумал, что это свойство довольно бесполезно, но, как оказалось, его можно применять для решения многих комбинаторных задач, в которых требуется вычислить формулу, задающуюся некоторым полиномом.↵

А теперь поставим общую задачу: пусть нам известны несколько значений, вычисленных в $k$ подряд идущих целых точках $x_i = i$ ($0 \leq i < k$), некоторой функции $f(x)$, являющейся многочленом. Требуется, зная эти точки, уметь находить значение $f(x)$ в произвольной точке $x'$. При этом гарантируется, что для старшей степени $d$ данного многочлена выполняется неравенство $d < k$ (сама эта задача где-то может называться как "задача об интерполяционном многочлене").↵

**Примечание**. В поставленной задаче **важно**, чтобы данная функция $f(x)$ являлась многочленом. Так, для функций, содержащих в себе операции, отличные от сложения, умножения и возведения в константную натуральную степень, приведённый ниже метод работать **не будет**.↵

Кажется, что данная задача решается только решением системы из $k$ линейных уравнений, где каждое уравнение &mdash; это многочлен степени $k-1$ с неизвестными коэффициентами. На самом деле существует множество различных способов решать эту задачу и в этом посте мы разберём один из них.↵

Так как же всё-таки находить многочлен, задающий $f(x)$, довольно просто и быстро и причём тут биномиальные коэффициенты?↵

**Примечание**. Биномиальный коэффициент $\binom{n}{k}$ &mdash; это то же самое, что и число сочетаний $C_n^k$. Они так называются, потому что они появляются при раскрытии бинома Ньютона.↵

Так получается, что любой многочлен $P(x)$ степени $d$ представляется в виде линейной комбинации $k+1$ биномиальных коэффициентов $\binom{x}{i}$, причём $0 \leq i \leq k$ и $k \geq d$.↵

Другими словами, любой полином $P(x)$ можно представить в виде суммы $\sum_{i=0}^{k} a_i \binom{x}{i}$, где $a_i$ &mdash; некоторые коэффициенты, причём в единственном виде (а для целых многочленов также верно, что $a_i$ &mdash; целые числа). Доказательство этого факта предоставляется читателю.↵

Осталось понять только то, чему равны коэффициенты $a_i$ в этой линейной комбинации. Оказывается, что $a_i = p_i(0)$, где $p_i$ &mdash; массив $i$-й разности исходного массива (что ещё иногда также называют "$i$-й производной на массиве"). Для случая $i = 0$ полагаем $p_0(j) = f(x_j)$, а при $i > 0$ $p_i$ определяем как $p_i(j) = p_{i-1}(j+1) - p_{i-1}(j)$ (если вы знаете, как доказывать этот факт, и не против поделиться, то напишите доказательство в комментариях).↵

<spoiler summary="Спойлер">↵
Для значений $p_i(0)$ также справедлива довольно легко выводящаяся явная формула (идея очень сильно напоминает разложение бинома Ньютона):  $p_i(0) = \sum_{j=0}^i (-1)^{i-j} \binom{i}{j} f(x_j)$. Её также можно использовать для достижения асимптотики, которую мы получим позже.↵
</spoiler>↵

Полученные массивы разности $p$ во многом схожи с привычными нам из математического анализа производными (а само разложение многочлена в линейную комбинацию напоминает ряд Тейлора). Приведу несколько их свойств:↵

- Для некоторого $j$ массив $p_j$ полностью состоит из нулей (или пуст) тогда и только тогда, когда искомый многочлен имеет степень $d < j$ (аналогично с производной функции, которая просто равна нулю).↵

- Каждый следующий массив разностей $p_j$ содержит на один элемент меньше, чем предыдущий $p_{j-1}$.↵

Из первого свойства следует, что если в задаче требуется найти только лишь степень многочлена, то достаточно вычислить массивы разностей и найти массив $p_i$ с минимальным $i$, состоящий из всех нулей (или пустой). Тогда степень искомого многочлена будет равна $i-1$ (конечно, для найденного многочлена $P(x)$ не всегда будет верно то, что он будет равен искомому многочлену $f(x)$, так как последний может иметь большую степень. Тогда нужно просто увеличить $k$, но обычно для многих комбинаторных задач с полиномиальными формулами старшая степень не превосходит $4$). ↵

<spoiler summary="Пример">↵
Пусть известны 5 элементов некоторого массива $a: [1;\;2;\;4;\;8;\;15]$. Тогда массивы разностей для него будут выглядеть так:↵

$$↵
\begin{matrix}↵
p_0: & \textbf{1} && 2 && 4 && 8 && 15\\↵
p_1: && \textbf{1} && 2 && 4 && 7\\↵
p_2: &&& \textbf{1} && 2 && 3\\↵
p_3: &&&& \textbf{1} && 1\\↵
p_4: &&&&& \textbf{0}\\↵
\end{matrix}↵
$$↵

Жирным шрифтом выделены элементы $p_i(0)$. Из таблицы разностей легко видеть, что $p_4$ состоит из нулей, а $p_3$ &mdash; нет $\implies$ многочлен, задающий массив $a$, имеет старшую степень $d = 3$. ↵

Сам многочлен имеет вид $P(x) = 1 \cdot \binom{x}{0} + 1 \cdot \binom{x}{1} + 1 \cdot \binom{x}{2} + 1 \cdot \binom{x}{3} + 0 \cdot \binom{x}{4} = \frac{1}{6}x^3 + \frac{5}{6}x + 1$. Нетрудно видеть, что $P(i) = a_i$, а старшая степень действительно равна $3$.↵
</spoiler>↵

Зная всю теорию, приведённую выше, уже можно написать алгоритм, работающий за $\mathcal{O}(k^2)$ времени и памяти, где $k$ &mdash; количество точек, для которых известно значение функции $f$. Имея величины $f(x_i)$, вычислим все значения $p_i(0)$ за $\mathcal{O}(k^2)$ времени и памяти (так как значения массива $i$-х разностей зависят только от массива $(i-1)$-х разностей, вообще можно хранить всего лишь два слоя, тогда потребление памяти сократится до $\mathcal{O}(k)$, но это не обязательно, ведь построение и так работает за $\mathcal{O}(k^2)$ времени). ↵

Если нам не нужно выводить формулу в явном виде для функции $f(x)$, то этого уже достаточно, чтобы решить задачу. Для любого запроса "узнать значение в какой-либо точке" можно просто посчитать значение суммы $\sum_{i=0}^{k} p_i(0) \binom{x}{i}$ за $\mathcal{O}(k^2\log \operatorname{MOD})$ (так как вычисления скорее всего производятся по модулю, то при подсчете $\binom{x}{k}$ вместо деления придётся $k$ раз умножать на обратный элемент за $\mathcal{O}(\log \operatorname{MOD})$). Но на самом деле можно быстрее!↵

Заметим, что при вычислении формулы $\sum_{i=0}^{k} p_i(0) \binom{x}{i}$ в её соседних членах нижняя часть соответствующих биноминальных коэффициентов отличается на 1. А так как для соседних биноминальных коэффициентов выполняется свойство $\binom{n}{k} = \frac{n-k+1}{k} \binom{n}{k-1}$, то подсчёт значения формулы можно ускорить до $\mathcal{O}(k \log \operatorname{MOD})$, причём код становится проще и короче, так как не нужно реализовывать функцию для вычисления $\binom{n}{k}$.↵

Таким образом, алгоритм работает за время $\mathcal{O}(k^2 + q k \log \operatorname{MOD})$, где $q$ &mdash; количество запросов на вычисление значения функции $f$ для какого-либо аргумента.↵

<spoiler summary="Код.">↵
~~~~~↵
#include <bits/stdc++.h>↵

const int MOD = 1e9 + 7;↵
//Вычисление
**UPDATE**. Алгоритм можно ещё больше упростить и ускорить, если заметить, что при подсчёте значения в какой-либо точке мы используем только обратногоые элемента по модулю↵
int inv(int a) {↵
int res = 1, p = MOD - 2;↵
while (p) {↵
if (p & 1)↵
res = 1LL * res * a % MOD;↵
a = 1LL * a * a % MOD;↵
p >>= 1;↵
}↵
return res;↵
}
ы для чисел $1,\,2,\,\dots,\,k$. Следовательно, можно в самом начале преподсчитать их значения, а потом уже использовать результаты. Таким образом, к времени предпосчёта добавляется член $\mathcal{O}(k)$, который не меняет асимптотику, а ответ на запрос можно теперь совершать за время $\mathcal{O}(k)$ (если предпосчитывать обратные элементы [за линейное время](https://e-maxx.ru/algo/reverse_element#4)). Общая асимптотика принимает вид $\mathcal{O}(k^2 + q k)$.↵

<spoiler summary="Код.">↵
~~~~~↵
#include <bits/stdc++.h>↵

const int MOD = 1e9 + 7;

struct PolynomialBuilder {↵
std::vector<std::vector<int>> p; ↵
std::vector<int> inv;↵
//Вычисление в конструкторе массивов разностей p↵
PolynomialBuilder(const std::vector<int> &x) {↵
p.push_back(x);↵
while ((int)p.back().size() > 1) {↵
std::vector<int> temp(p.back().size() - 1);↵
for (int i = 0; i < (int)temp.size(); i++)↵
temp[i] = (p.back()[i + 1] + MOD - p.back()[i]) % MOD;↵
p.push_back(temp);↵
}↵
//Предпосчёт обратных элементов↵
inv.resize(p.size() + 1);↵
inv[1] = 1;↵
for (int i = 2; i < (int)inv.size(); i++)↵
inv[i] = (MOD - 1LL * (MOD / i) * inv[MOD % i] % MOD) % MOD;↵
}↵
//Вычисление значения в точке для значений |x| ~ 1e18↵
int evaluate_at(long long x) { ↵
x = (x % MOD + MOD) % MOD;↵
int res = 0, binom = 1;↵
for (int i = 0; i < (int)p.size(); i++) {↵
//Если n < k и n >= 0, то Cnk(n, k) = 0↵
if (x >= 0 && x < i)↵
break;
 
res = (res + 1LL * p[i][0] * binom % MOD + MOD) % MOD;↵
binom = 1LL * binom * (x - i) % MOD * inv
([i + 1)] % MOD;↵
}↵
return res;↵
}↵
};↵

int main() {↵
PolynomialBuilder f({0, 1, 4});  //Задаем точками полином x^2↵
std::cout << f.evaluate_at(3) << " ";  //Вывод: 9↵
std::cout << f.evaluate_at(-5);  //Вывод: 25↵
}↵

~~~~~↵
</spoiler>↵

Но у приведённого выше алгоритма есть одна проблема: он работает за приведённую асимптотику, если мы уже **преподсчитали** значения $f(x_i)$ (то есть если значения в точках $x_i$ каким-либо образом зависят от входных данных, то предпосчитать вне программы уже не получится). ↵

Эта проблема решается внедрением предпосчёта в программу. Если вычисление значений функции $f$ в точках $x_i$ занимает $\mathcal{O}(m)$ времени, то общая асимптотика примет вид $\mathcal{O}(mk + k^2 + q k \log \operatorname{MOD})$ времени. Здесь становится видно, что для решения комбинаторной задачи достаточно запустить какое-либо простое и медленное решение $k$ раз, которое посчитает значение функции в начальных точках. И очень часто даже экспоненциальные алгоритмы подходят для предпосчета, так как в итоговом многочлене степень обычно достаточно маленькая.↵

Если же в задаче требуется вывести явную формулу, вычисляющую функцию, то это уже будет не так просто реализовать, ведь придётся использовать реализацию рациональных чисел для коэффициентов (Fraction в python или самописный для малых $k$), хранить многочлены в какой-либо структуре данных (можно хранить коэффициенты в std::vector из C++), а также уметь их складывать и умножать. Все это реально реализовать и заставить работать за асимптотику $\mathcal{O}(mk + k^3 + qk)$ с получением ответа на запрос "узнать значение в какой-либо точке $n$" за $\mathcal{O}(k)$ (но уже для $k > 20$ придётся использовать длинную арифметику, что влечёт за собой серьезное увеличение времени работы алгоритма, так как знаменатели коэффициентов многочлена могут стать порядка $k! \approx 10^{18}$). ↵

Таким образом, метод, приведённый в этом посте, довольно быстро решает задачу об интерполяционном многочлене. Так, если многочлен не нужно выводить в явном виде, то он работает быстрее, чем работал бы метод Гаусса, не говоря уже о простоте реализации.↵

Список задач, на которых можно опробовать этот метод (если вы знаете какие-либо ещё задачи, то можете их оставлять в комментариях, я их буду добавлять сюда):↵

- [Задача G "Пазл" с финала ВКОШП 2017](https://mirror.codeforces.com/gym/101636)↵

- [Задача F "Fygon" с NWERC ICPC 2015](http://neerc.itmo.ru/archive/2015.html)↵

- [Продолжи последовательность (Timus №2125)](https://timus.online/problem.aspx?space=1&num=2125) &mdash; приведённый в посте алгоритм не зайдёт, но, возможно, существует какая-то модификация алгоритма, работающая за $\mathcal{O}(k \log k \log \operatorname{MOD})$ (если есть идеи по модификации, буду рад их услышать в комментариях).

History

 
 
 
 
Revisions
 
 
  Rev. Lang. By When Δ Comment
ru9 Russian unreal.eugene 2019-10-23 12:48:21 24
ru8 Russian unreal.eugene 2019-10-23 12:40:39 280
ru7 Russian unreal.eugene 2019-10-23 12:38:15 1036
ru6 Russian unreal.eugene 2019-10-23 10:17:48 71 Мелкая правка: 'к:\n\n$$\np0:124815p1:1247p2:123p3:11p4:0\begin{mat' -> 'к:\n\n$$\n\begin{mat'
ru5 Russian unreal.eugene 2019-10-23 10:02:55 3
ru4 Russian unreal.eugene 2019-10-23 09:54:58 9 (опубликовано)
ru3 Russian unreal.eugene 2019-10-23 09:51:29 8
ru2 Russian unreal.eugene 2019-10-23 09:49:24 233
ru1 Russian unreal.eugene 2019-10-23 01:18:45 11225 Первая редакция (сохранено в черновиках)