Привет всем! ↵
↵
Года два назад я решал задачу, в которой для сдачи требовалось вывести формулу, сворачивая суммы, ряды и делая кучу неприятных дел, на которые уходит много сил и энергии на контесте. Я тогда заинтересовался, а существует ли более простой способ вывести эту формулу, не делая почти никаких математических преобразований?↵
↵
[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$ линейных уравнений, где каждое уравнение — это многочлен степени $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)$ (если вы знаете, как доказывать этот факт, и не против поделиться, то напишите доказательство в комментариях).↵
↵
<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$ — нет $\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$ — количество точек, для которых известно значение функции $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$ для какого-либо аргумента.↵
↵
<spoiler summary="Код.">↵
~~~~~↵
#include <bits/stdc++.h>↵
↵
const int MOD = 1e9 + 7;↵
//Вычисление обратного элемента по модулю↵
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;↵
}↵
struct PolynomialBuilder {↵
std::vector<std::vector<int>> p; ↵
//Вычисление в конструкторе массивов разностей 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);↵
}↵
}↵
//Вычисление значения в точке для значений |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++) {↵
if (x >= 0 && x < i)↵
break;//Если 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://codeforces.net/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) — приведённый в посте алгоритм не зайдёт, но, возможно, существует какая-то модификация алгоритма, работающая за $\mathcal{O}(k \log k \log \operatorname{MOD})$ (если есть идеи по модификации, буду рад их услышать в комментариях).
↵
Года два назад я решал задачу, в которой для сдачи требовалось вывести формулу, сворачивая суммы, ряды и делая кучу неприятных дел, на которые уходит много сил и энергии на контесте. Я тогда заинтересовался, а существует ли более простой способ вывести эту формулу, не делая почти никаких математических преобразований?↵
↵
[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$ линейных уравнений, где каждое уравнение — это многочлен степени $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)$ (если вы знаете, как доказывать этот факт, и не против поделиться, то напишите доказательство в комментариях).↵
↵
<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$ — нет $\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$ — количество точек, для которых известно значение функции $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$ для какого-либо аргумента.↵
↵
<spoiler summary="Код.">↵
~~~~~↵
#include <bits/stdc++.h>↵
↵
const int MOD = 1e9 + 7;↵
//Вычисление обратного элемента по модулю↵
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;↵
}↵
struct PolynomialBuilder {↵
std::vector<std::vector<int>> p; ↵
//Вычисление в конструкторе массивов разностей 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);↵
}↵
}↵
//Вычисление значения в точке для значений |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++) {↵
break;
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://codeforces.net/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) — приведённый в посте алгоритм не зайдёт, но, возможно, существует какая-то модификация алгоритма, работающая за $\mathcal{O}(k \log k \log \operatorname{MOD})$ (если есть идеи по модификации, буду рад их услышать в комментариях).
