Привет всем!
Года два назад я решал задачу, в которой для сдачи требовалось вывести формулу, сворачивая суммы, ряды и делая кучу неприятных дел, на которые уходит много сил и энергии на контесте. Я тогда заинтересовался, а существует ли более простой способ вывести эту формулу, не делая почти никаких математических преобразований? Задолго до этого я наткнулся на одно свойство биномиальных коэффициентов в Википедии. Я естественно подумал, что это свойство довольно бесполезно, но, как оказалось, его можно применять для решения многих комбинаторных задач, в которых требуется вычислить формулу, задающуюся некоторым полиномом.
А теперь поставим общую задачу: пусть нам известны несколько значений, вычисленных в $$$k$$$ подряд идущих целых точках $$$x_i = i$$$ ($$$0 \leq i < k$$$), некоторой функции $$$f(x)$$$, являющейся многочленом. Требуется, зная эти точки, уметь находить значение $$$f(x)$$$ в произвольной точке $$$x'$$$. При этом гарантируется, что для старшей степени $$$d$$$ данного многочлена выполняется неравенство $$$d < k$$$ (сама эта задача где-то может называться как "задача об интерполяционном многочлене").
Примечание. В поставленной задаче важно, чтобы данная функция $$$f(x)$$$ являлась многочленом. Так, для функций, содержащих в себе операции, отличные от сложения, умножения и возведения в константную натуральную степень, приведённый ниже метод работать не будет.
Кажется, что данная задача решается только решением системы из $$$k$$$ линейных уравнений, где каждое уравнение — это многочлен степени $$$k-1$$$ с неизвестными коэффициентами. На самом деле существует множество различных способов решать эту задачу и в этом посте мы разберём один из них.
Так как же всё-таки находить многочлен, задающий $$$f(x)$$$, довольно просто и быстро и причём тут биномиальные коэффициенты?
Примечание. Биномиальный коэффициент $$$\binom{n}{k}$$$ — это то же самое, что и число сочетаний $$$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$$$ — некоторые коэффициенты, причём в единственном виде (а для целых многочленов также верно, что $$$a_i$$$ — целые числа). Доказательство этого факта предоставляется читателю.
Осталось понять только то, чему равны коэффициенты $$$a_i$$$ в этой линейной комбинации. Оказывается, что $$$a_i = p_i(0)$$$, где $$$p_i$$$ — массив $$$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)$$$ (если вы знаете, как доказывать этот факт, и не против поделиться, то напишите доказательство в комментариях).
Полученные массивы разности $$$p$$$ во многом схожи с привычными нам из математического анализа производными (а само разложение многочлена в линейную комбинацию напоминает ряд Тейлора). Приведу несколько их свойств:
Для некоторого $$$j$$$ массив $$$p_j$$$ полностью состоит из нулей (или пуст) тогда и только тогда, когда искомый многочлен имеет степень $$$d < j$$$ (аналогично с производной функции, которая просто равна нулю).
Каждый следующий массив разностей $$$p_j$$$ содержит на один элемент меньше, чем предыдущий $$$p_{j-1}$$$.
Из первого свойства следует, что если в задаче требуется найти только лишь степень многочлена, то достаточно вычислить массивы разностей и найти массив $$$p_i$$$ с минимальным $$$i$$$, состоящий из всех нулей (или пустой). Тогда степень искомого многочлена будет равна $$$i-1$$$ (конечно, для найденного многочлена $$$P(x)$$$ не всегда будет верно то, что он будет равен искомому многочлену $$$f(x)$$$, так как последний может иметь большую степень. Тогда нужно просто увеличить $$$k$$$, но обычно для многих комбинаторных задач с полиномиальными формулами старшая степень не превосходит $$$4$$$).
Зная всю теорию, приведённую выше, уже можно написать алгоритм, работающий за $$$\mathcal{O}(k^2)$$$ времени и памяти, где $$$k$$$ — количество точек, для которых известно значение функции $$$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$$$ — количество запросов на вычисление значения функции $$$f$$$ для какого-либо аргумента.
Но у приведённого выше алгоритма есть одна проблема: он работает за приведённую асимптотику, если мы уже преподсчитали значения $$$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}$$$).
Таким образом, метод, приведённый в этом посте, довольно быстро решает задачу об интерполяционном многочлене. Так, если многочлен не нужно выводить в явном виде, то он работает быстрее, чем работал бы метод Гаусса, не говоря уже о простоте реализации.
Список задач, на которых можно опробовать этот метод (если вы знаете какие-либо ещё задачи, то можете их оставлять в комментариях, я их буду добавлять сюда):
Продолжи последовательность (Timus №2125) — приведённый в посте алгоритм не зайдёт, но, возможно, существует какая-то модификация алгоритма, работающая за $$$\mathcal{O}(k \log k \log \operatorname{MOD})$$$ (если есть идеи по модификации, буду рад их услышать в комментариях).
