Операції над многочленами та рядами

Задачі зі змагального програмування, особливо ті, що так чи інакше пов'язані з перелічуванням, часто розв'язуються зведенням задачі до обчислення чогось над многочленами та формальними степеневими рядами.

Сюди входять такі поняття, як множення многочленів, інтерполяція, а також складніші — наприклад, логарифми та експоненти многочленів. У цій статті подано короткий огляд таких операцій та поширених підходів до них.

Базові поняття та факти

У цьому розділі ми зосереджуємося радше на означеннях та «інтуїтивних» властивостях різних операцій над многочленами. Технічні деталі їхньої реалізації та складності будуть розглянуті в подальших розділах.

Множення многочленів

Означення

Многочлен від однієї змінної — це вираз вигляду $A(x) = a_0 + a_1 x + \dots + a_n x^n$.

Значення $a_0, \dots, a_n$ — це коефіцієнти многочлена, які зазвичай беруться з деякої множини чисел або числоподібних структур. У цій статті ми припускаємо, що коефіцієнти беруться з деякого поля, тобто операції додавання, віднімання, множення та ділення для них коректно визначені (окрім ділення на $0$) і вони загалом поводяться подібно до дійсних чисел.

Типовий приклад такого поля — поле остач за простим модулем $p$.

Для простоти ми опускатимемо слова від однієї змінної, оскільки це єдиний вид многочленів, який ми розглядаємо в цій статті. Ми також писатимемо $A$ замість $A(x)$ скрізь, де це можливо, що буде зрозуміло з контексту. Припускається, що або $a_n \neq 0$, або $A(x)=0$.

Означення

Добутком двох многочленів є вираз, отриманий розкриттям дужок як арифметичного виразу:

$$A(x) B(x) = \left(\sum\limits_{i=0}^n a_i x^i \right)\left(\sum\limits_{j=0}^m b_j x^j\right) = \sum\limits_{i,j} a_i b_j x^{i+j} = \sum\limits_{k=0}^{n+m} c_k x^k = C(x).$$

Послідовність $c_0, c_1, \dots, c_{n+m}$ коефіцієнтів $C(x)$ називається згорткою послідовностей $a_0, \dots, a_n$ та $b_0, \dots, b_m$.

Означення

Степінь многочлена $A$ з $a_n \neq 0$ визначається як $\deg A = n$.

Для узгодженості степінь $A(x) = 0$ визначається як $\deg A = -\infty$.

У такому означенні $\deg AB = \deg A + \deg B$ для будь-яких многочленів $A$ та $B$.

Згортки лежать в основі розв'язання багатьох перелічувальних задач.

інформація

У вас є $n$ об'єктів першого типу та $m$ об'єктів другого типу.

Об'єкти першого типу мають значення $a_1, \dots, a_n$, а об'єкти другого типу — значення $b_1, \dots, b_m$.

Ви вибираєте один об'єкт першого типу та один об'єкт другого типу. Скількома способами можна отримати сумарне значення $k$?

Розв'язок

Розглянемо добуток $(x^{a_1} + \dots + x^{a_n})(x^{b_1} + \dots + x^{b_m})$. Якщо розкрити дужки, кожен одночлен відповідатиме парі $(a_i, b_j)$ і даватиме внесок у коефіцієнт біля $x^{a_i+b_j}$. Іншими словами, відповіддю є коефіцієнт біля $x^k$ у добутку.

інформація

Ви кидаєте $6$-гранний кубик $n$ разів і підсумовуєте результати всіх кидків. Яка ймовірність отримати суму $k$?

Розв'язок

Відповідь — це кількість результатів, що мають суму $k$, поділена на загальну кількість результатів, яка дорівнює $6^n$.

Яка кількість результатів, що мають суму $k$? Для $n=1$ її можна подати многочленом $A(x) = x^1+x^2+\dots+x^6$.

Для $n=2$, використовуючи той самий підхід, що й у прикладі вище, ми робимо висновок, що вона подається многочленом $(x^1+x^2+\dots+x^6)^2$.

Зважаючи на це, відповіддю на задачу є $k$-й коефіцієнт многочлена $(x^1+x^2+\dots+x^6)^n$, поділений на $6^n$.

Коефіцієнт біля $x^k$ у многочлені $A(x)$ коротко позначається як $[x^k]A$.

Формальний степеневий ряд

Означення

Формальний степеневий ряд — це нескінченна сума $A(x) = a_0 + a_1 x + a_2 x^2 + \dots$, що розглядається безвідносно до її властивостей збіжності.

Іншими словами, коли ми розглядаємо, наприклад, суму $1+\frac{1}{2}+\frac{1}{4}+\frac{1}{8}+\dots=2$, ми маємо на увазі, що вона збігається до $2$, коли кількість доданків прямує до нескінченності. Однак формальні ряди розглядаються лише в термінах послідовностей, які їх утворюють.

Означення

Добуток формальних степеневих рядів $A(x)$ та $B(x)$ також визначається розкриттям дужок як арифметичного виразу:

$$A(x) B(x) = \left(\sum\limits_{i=0}^\infty a_i x^i \right)\left(\sum\limits_{j=0}^\infty b_j x^j\right) = \sum\limits_{i,j} a_i b_j x^{i+j} = \sum\limits_{k=0}^{\infty} c_k x^k = C(x),$$

де коефіцієнти $c_0, c_1, \dots$ визначаються як скінченні суми

$$c_k = \sum\limits_{i=0}^k a_i b_{k-i}.$$

Послідовність $c_0, c_1, \dots$ також називається згорткою послідовностей $a_0, a_1, \dots$ та $b_0, b_1, \dots$, узагальнюючи це поняття на нескінченні послідовності.

Отже, многочлени можна розглядати як формальні степеневі ряди, але зі скінченною кількістю коефіцієнтів.

Формальні степеневі ряди відіграють вирішальну роль у перелічувальній комбінаториці, де їх вивчають як твірні функції для різних послідовностей. Детальне пояснення твірних функцій та інтуїції за ними, на жаль, виходить за межі цієї статті, тому допитливого читача відсилаємо, наприклад, сюди за деталями про їхній комбінаторний зміст.

Однак ми дуже стисло зазначимо, що якщо $A(x)$ та $B(x)$ — це твірні функції для послідовностей, що перелічують деякі об'єкти за кількістю «атомів» у них (наприклад, дерева за кількістю вершин), то добуток $A(x) B(x)$ перелічує об'єкти, які можна описати як пари об'єктів типів $A$ та $B$, перелічуючи їх за загальною кількістю «атомів» у парі.

інформація

Нехай $A(x) = \sum\limits_{i=0}^\infty 2^i x^i$ перелічує набори каменів, кожен камінь у яких пофарбований в один з $2$ кольорів (тож існує $2^i$ таких наборів розміру $i$), а $B(x) = \sum\limits_{j=0}^{\infty} 3^j x^j$ перелічує набори каменів, кожен камінь у яких пофарбований в один з $3$ кольорів. Тоді $C(x) = A(x) B(x) = \sum\limits_{k=0}^\infty c_k x^k$ перелічуватиме об'єкти, які можна описати як «два набори каменів, перший набір лише з каменів типу $A$, другий набір лише з каменів типу $B$, із загальною кількістю каменів $k$» для $c_k$.

Подібним чином є інтуїтивний зміст і в деяких інших функцій над формальними степеневими рядами.

Ділення многочленів стовпчиком

Подібно до цілих чисел, для многочленів можна визначити ділення стовпчиком.

Означення

Для будь-яких многочленів $A$ та $B \neq 0$ можна подати $A$ у вигляді

$$A = D \cdot B + R,~ \deg R < \deg B,$$

де $R$ називається остачею від ділення $A$ на $B$, а $D$ називається часткою.

Позначаючи $\deg A = n$ та $\deg B = m$, наївний спосіб зробити це — застосувати ділення стовпчиком, під час якого ми множимо $B$ на одночлен $\frac{a_n}{b_m} x^{n - m}$ і віднімаємо результат від $A$, доки степінь $A$ не стане меншим за степінь $B$. Те, що залишиться від $A$ наприкінці, буде остачею (звідси й назва), а многочлени, на які ми множили $B$ у процесі, підсумовані разом, утворюють частку.

Означення

Якщо $A$ та $B$ мають однакову остачу за модулем $C$, то кажуть, що вони еквівалентні за модулем $C$, що позначається як

$$A \equiv B \pmod{C}.$$

Ділення многочленів стовпчиком корисне завдяки своїм численним важливим властивостям:

• $A$ кратне $B$ тоді й лише тоді, коли $A \equiv 0 \pmod B$.

• Звідси випливає, що $A \equiv B \pmod C$ тоді й лише тоді, коли $A-B$ кратне $C$.

• Зокрема, $A \equiv B \pmod{C \cdot D}$ тягне за собою $A \equiv B \pmod{C}$.

• Для будь-якого лінійного многочлена $x-r$ виконується $A(x) \equiv A(r) \pmod{x-r}$.

• Звідси випливає, що $A$ кратне $x-r$ тоді й лише тоді, коли $A(r)=0$.

• Для модуля $x^k$ виконується $A \equiv a_0 + a_1 x + \dots + a_{k-1} x^{k-1} \pmod{x^k}$.

Зауважимо, що ділення стовпчиком не можна коректно визначити для формальних степеневих рядів. Натомість для будь-якого $A(x)$ такого, що $a_0 \neq 0$, можна визначити обернений формальний степеневий ряд $A^{-1}(x)$ такий, що $A(x) A^{-1}(x) = 1$. Цей факт, своєю чергою, можна використати для обчислення результату ділення многочленів стовпчиком.

Базова реалізація

Тут ви можете знайти базову реалізацію алгебри многочленів.

Вона підтримує всі тривіальні операції та деякі інші корисні методи. Головний клас — poly<T> для многочленів з коефіцієнтами типу T.

Підтримуються всі арифметичні операції +, -, *, % та /, де % та / позначають остачу та частку в евклідовому діленні.

Також є клас modular<m> для виконання арифметичних операцій над остачами за простим модулем m.

Інші корисні функції:

• deriv(): обчислює похідну $P'(x)$ від $P(x)$.
• integr(): обчислює невизначений інтеграл $Q(x) = \int P(x)$ від $P(x)$ такий, що $Q(0)=0$.
• inv(size_t n): обчислює перші $n$ коефіцієнтів $P^{-1}(x)$ за $O(n \log n)$.
• log(size_t n): обчислює перші $n$ коефіцієнтів $\ln P(x)$ за $O(n \log n)$.
• exp(size_t n): обчислює перші $n$ коефіцієнтів $\exp P(x)$ за $O(n \log n)$.
• pow(size_t k, size_t n): обчислює перші $n$ коефіцієнтів $P^{k}(x)$ за $O(n \log nk)$.
• deg(): повертає степінь $P(x)$.
• lead(): повертає коефіцієнт біля $x^{\deg P(x)}$.
• resultant(poly<T> a, poly<T> b): обчислює результант $a$ та $b$ за $O(|a| \cdot |b|)$.
• bpow(T x, size_t n): обчислює $x^n$.
• bpow(T x, size_t n, T m): обчислює $x^n \pmod{m}$.
• chirpz(T z, size_t n): обчислює $P(1), P(z), P(z^2), \dots, P(z^{n-1})$ за $O(n \log n)$.
• vector<T> eval(vector<T> x): обчислює $P(x_1), \dots, P(x_n)$ за $O(n \log^2 n)$.
• poly<T> inter(vector<T> x, vector<T> y): інтерполює многочлен за набором пар $P(x_i) = y_i$ за $O(n \log^2 n)$.
• І ще деякі — не соромтеся досліджувати код!

Арифметика

Множення

Найголовніша операція — це множення двох многочленів. Тобто, маючи многочлени $A$ та $B$:

$$A = a_0 + a_1 x + \dots + a_n x^n$$ $$B = b_0 + b_1 x + \dots + b_m x^m$$

потрібно обчислити многочлен $C = A \cdot B$, який визначається як

$$\boxed{C = \sum\limits_{i=0}^n \sum\limits_{j=0}^m a_i b_j x^{i+j}} = c_0 + c_1 x + \dots + c_{n+m} x^{n+m}.$$

Його можна обчислити за $O(n \log n)$ за допомогою швидкого перетворення Фур'є, і майже всі методи тут використовуватимуть його як підпрограму.

Обернений ряд

Якщо $A(0) \neq 0$, то завжди існує нескінченний формальний степеневий ряд $A^{-1}(x) = q_0+q_1 x + q_2 x^2 + \dots$ такий, що $A^{-1} A = 1$. Часто буває корисно обчислити перші $k$ коефіцієнтів $A^{-1}$ (тобто обчислити його за модулем $x^k$). Є два основні способи зробити це.

«Розділяй і володарюй»

Цей алгоритм був згаданий у статті Шьонхаге і натхненний методом Греффе. Відомо, що для $B(x)=A(x)A(-x)$ виконується $B(x)=B(-x)$, тобто $B(x)$ — парний многочлен. Це означає, що він має ненульові коефіцієнти лише при парних номерах і може бути поданий як $B(x)=T(x^2)$. Таким чином, ми можемо зробити такий перехід:

$$A^{-1}(x) \equiv \frac{1}{A(x)} \equiv \frac{A(-x)}{A(x)A(-x)} \equiv \frac{A(-x)}{T(x^2)} \pmod{x^k}$$

Зауважимо, що $T(x)$ можна обчислити за допомогою одного множення, після чого нас цікавить лише перша половина коефіцієнтів його оберненого ряду. Це фактично зводить початкову задачу обчислення $A^{-1} \pmod{x^k}$ до обчислення $T^{-1} \pmod{x^{\lceil k / 2 \rceil}}$.

Складність цього методу можна оцінити як

$$T(n) = T(n/2) + O(n \log n) = O(n \log n).$$

Алгоритм Зівекінга–Кунга

Загальний процес, описаний тут, відомий як підняття Гензеля, оскільки він випливає з леми Гензеля. Ми розглянемо його детальніше далі, а поки що зосередимося на спеціальному розв'язку. Частина «підняття» тут означає, що ми починаємо з наближення $B_0=q_0=a_0^{-1}$, яке є $A^{-1} \pmod x$, а потім ітеративно піднімаємося від $\bmod x^a$ до $\bmod x^{2a}$.

Нехай $B_k \equiv A^{-1} \pmod{x^a}$. Наступне наближення має задовольняти рівняння $A B_{k+1} \equiv 1 \pmod{x^{2a}}$ і може бути подане як $B_{k+1} = B_k + x^a C$. Звідси випливає рівняння

$$A(B_k + x^{a}C) \equiv 1 \pmod{x^{2a}}.$$

Нехай $A B_k \equiv 1 + x^a D \pmod{x^{2a}}$, тоді рівняння вище тягне за собою

$$x^a(D+AC) \equiv 0 \pmod{x^{2a}} \implies D \equiv -AC \pmod{x^a} \implies C \equiv -B_k D \pmod{x^a}.$$

Звідси можна отримати остаточну формулу, яка має вигляд

$$x^a C \equiv -B_k x^a D \equiv B_k(1-AB_k) \pmod{x^{2a}} \implies \boxed{B_{k+1} \equiv B_k(2-AB_k) \pmod{x^{2a}}}$$

Отже, почавши з $B_0 \equiv a_0^{-1} \pmod x$, ми обчислимо послідовність $B_k$ таку, що $AB_k \equiv 1 \pmod{x^{2^k}}$, зі складністю

$$T(n) = T(n/2) + O(n \log n) = O(n \log n).$$

Алгоритм тут може здатися дещо складнішим за перший, але за ним стоїть дуже ґрунтовне й практичне обґрунтування, а також великий потенціал до узагальнення, якщо подивитися на нього з іншого боку, що буде пояснено далі.

Евклідове ділення

Розглянемо два многочлени $A(x)$ та $B(x)$ степенів $n$ та $m$. Як було сказано раніше, $A(x)$ можна переписати як

$$A(x) = B(x) D(x) + R(x), \deg R < \deg B.$$

Нехай $n \geq m$, це означало б, що $\deg D = n - m$ і старші $n-m+1$ коефіцієнтів $A$ не впливають на $R$. Це означає, що $D(x)$ можна відновити зі старших $n-m+1$ коефіцієнтів $A(x)$ та $B(x)$, якщо розглянути це як систему рівнянь.

Систему лінійних рівнянь, про яку йдеться, можна записати в такому вигляді:

$$\begin{bmatrix} a_n \\ \vdots \\ a_{m+1} \\ a_{m} \end{bmatrix} = \begin{bmatrix} b_m & \dots & 0 & 0 \\ \vdots & \ddots & \vdots & \vdots \\ \dots & \dots & b_m & 0 \\ \dots & \dots & b_{m-1} & b_m \end{bmatrix} \begin{bmatrix}d_{n-m} \\ \vdots \\ d_1 \\ d_0\end{bmatrix}$$

З цього вигляду можна зробити висновок, що з введенням обернених многочленів

$$A^R(x) = x^nA(x^{-1})= a_n + a_{n-1} x + \dots + a_0 x^n$$ $$B^R(x) = x^m B(x^{-1}) = b_m + b_{m-1} x + \dots + b_0 x^m$$ $$D^R(x) = x^{n-m}D(x^{-1}) = d_{n-m} + d_{n-m-1} x + \dots + d_0 x^{n-m}$$

систему можна переписати як

$$A^R(x) \equiv B^R(x) D^R(x) \pmod{x^{n-m+1}}.$$

Звідси можна однозначно відновити всі коефіцієнти $D(x)$:

$$\boxed{D^R(x) \equiv A^R(x) (B^R(x))^{-1} \pmod{x^{n-m+1}}}$$

А з цього, своєю чергою, можна відновити $R(x)$ як $R(x) = A(x) - B(x)D(x)$.

Зауважимо, що матриця вище — це так звана трикутна матриця Тепліца, і, як ми бачимо, розв'язання системи лінійних рівнянь з довільною матрицею Тепліца насправді еквівалентне оберненню многочлена. Більше того, обернена до неї матриця також буде трикутною матрицею Тепліца, а її елементи, у термінах, що використовувалися вище, є коефіцієнтами $(B^R(x))^{-1} \pmod{x^{n-m+1}}$.

Обчислення функцій від многочлена

Метод Ньютона

Узагальнимо алгоритм Зівекінга–Кунга. Розглянемо рівняння $F(P) = 0$, де $P(x)$ має бути многочленом, а $F(x)$ — деяка функція зі значеннями-многочленами, визначена як

$$F(x) = \sum\limits_{i=0}^\infty \alpha_i (x-\beta)^i,$$

де $\beta$ — деяка константа. Можна довести, що, ввівши нову формальну змінну $y$, ми можемо виразити $F(x)$ як

$$F(x) = F(y) + (x-y)F'(y) + (x-y)^2 G(x,y),$$

де $F'(x)$ — формальний степеневий ряд похідної, визначений як

$$F'(x) = \sum\limits_{i=0}^\infty (i+1)\alpha_{i+1}(x-\beta)^i,$$

а $G(x, y)$ — деякий формальний степеневий ряд від $x$ та $y$. З цим результатом ми можемо знайти розв'язок ітеративно.

Нехай $F(Q_k) \equiv 0 \pmod{x^{a}}$. Нам потрібно знайти $Q_{k+1} \equiv Q_k + x^a C \pmod{x^{2a}}$ такий, що $F(Q_{k+1}) \equiv 0 \pmod{x^{2a}}$.

Підставивши $x = Q_{k+1}$ та $y=Q_k$ у формулу вище, ми отримуємо

$$F(Q_{k+1}) \equiv F(Q_k) + (Q_{k+1} - Q_k) F'(Q_k) + (Q_{k+1} - Q_k)^2 G(x, y) \pmod x^{2a}.$$

Оскільки $Q_{k+1} - Q_k \equiv 0 \pmod{x^a}$, то також виконується $(Q_{k+1} - Q_k)^2 \equiv 0 \pmod{x^{2a}}$, тому

$$0 \equiv F(Q_{k+1}) \equiv F(Q_k) + (Q_{k+1} - Q_k) F'(Q_k) \pmod{x^{2a}}.$$

Остання формула дає нам значення $Q_{k+1}$:

$$\boxed{Q_{k+1} = Q_k - \dfrac{F(Q_k)}{F'(Q_k)} \pmod{x^{2a}}}$$

Отже, знаючи, як обертати многочлени та як обчислювати $F(Q_k)$, ми можемо знайти $n$ коефіцієнтів $P$ зі складністю

$$T(n) = T(n/2) + f(n),$$

де $f(n)$ — час, потрібний для обчислення $F(Q_k)$ та $F'(Q_k)^{-1}$, який зазвичай становить $O(n \log n)$.

Ітеративне правило вище відоме в чисельному аналізі як метод Ньютона.

Лема Гензеля

Як було згадано раніше, формально й у загальному вигляді цей результат відомий як лема Гензеля, і її насправді можна використовувати в ще ширшому сенсі, коли ми працюємо з послідовністю вкладених кілець. У цьому конкретному випадку ми працювали з послідовністю остач многочленів за модулями $x$, $x^2$, $x^3$ і так далі.

Інший приклад, де підняття Гензеля може бути корисним, — це так звані $p$-адичні числа, де ми фактично працюємо з послідовністю остач цілих чисел за модулями $p$, $p^2$, $p^3$ і так далі. Наприклад, метод Ньютона можна використати, щоб знайти всі можливі автоморфні числа (числа, які закінчуються самі на себе при піднесенні до квадрата) із заданою основою системи числення. Цю задачу залишаємо читачеві як вправу. Ви можете розглянути цю задачу, щоб перевірити, чи працює ваш розв'язок для чисел з основою $10$.

Логарифм

Для функції $\ln P(x)$ відомо, що:

$$\boxed{(\ln P(x))' = \dfrac{P'(x)}{P(x)}}$$

Таким чином, ми можемо обчислити $n$ коефіцієнтів $\ln P(x)$ за $O(n \log n)$.

Обернений ряд

Виявляється, ми можемо отримати формулу для $A^{-1}$ за допомогою методу Ньютона. Для цього ми беремо рівняння $A=Q^{-1}$, отже:

$$F(Q) = Q^{-1} - A$$ $$F'(Q) = -Q^{-2}$$ $$\boxed{Q_{k+1} \equiv Q_k(2-AQ_k) \pmod{x^{2^{k+1}}}}$$

Експонента

Навчимося обчислювати $e^{P(x)}=Q(x)$. Має виконуватися $\ln Q = P$, отже:

$$F(Q) = \ln Q - P$$ $$F'(Q) = Q^{-1}$$ $$\boxed{Q_{k+1} \equiv Q_k(1 + P - \ln Q_k) \pmod{x^{2^{k+1}}}}$$

$k$-й степінь

Тепер нам потрібно обчислити $P^k(x)=Q$. Це можна зробити за допомогою такої формули:

$$Q = \exp\left[k \ln P(x)\right]$$

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

Щоб його знайти, потрібно обчислити логарифм або експоненту вільного коефіцієнта многочлена.

Але єдиний розумний спосіб це зробити — якщо $P(0)=1$ для $Q = \ln P$, тоді $Q(0)=0$, і якщо $P(0)=0$ для $Q = e^P$, тоді $Q(0)=1$.

Отже, формулу вище можна використовувати лише якщо $P(0) = 1$. Інакше, якщо $P(x) = \alpha x^t T(x)$, де $T(0)=1$, можна записати, що:

$$\boxed{P^k(x) = \alpha^kx^{kt} \exp[k \ln T(x)]}$$

Зауважте, що ви також можете обчислити деякий корінь $k$-го степеня з многочлена, якщо можете обчислити $\sqrt[k]{\alpha}$, наприклад для $\alpha=1$.

Обчислення в точках та інтерполяція

Chirp-z перетворення

Для окремого випадку, коли вам потрібно обчислити многочлен у точках $x_r = z^{2r}$, можна зробити таке:

$$A(z^{2r}) = \sum\limits_{k=0}^n a_k z^{2kr}$$

Підставимо $2kr = r^2+k^2-(r-k)^2$. Тоді ця сума переписується як:

$$\boxed{A(z^{2r}) = z^{r^2}\sum\limits_{k=0}^n (a_k z^{k^2}) z^{-(r-k)^2}}$$

Що, з точністю до множника $z^{r^2}$, дорівнює згортці послідовностей $u_k = a_k z^{k^2}$ та $v_k = z^{-k^2}$.

Зауважте, що тут $u_k$ має індекси від $0$ до $n$, а $v_k$ має індекси від $-n$ до $m$, де $m$ — максимальний степінь $z$, який вам потрібен.

Тепер, якщо вам потрібно обчислити многочлен у точках $x_r = z^{2r+1}$, ви можете звести це до попередньої задачі за допомогою перетворення $a_k \to a_k z^k$.

Це дає нам алгоритм за $O(n \log n)$, коли вам потрібно обчислити значення в степенях $z$, тож ви можете обчислити DFT для розмірів, що не є степенями двійки.

Інше спостереження полягає в тому, що $kr = \binom{k+r}{2} - \binom{k}{2} - \binom{r}{2}$. Тоді маємо

$$\boxed{A(z^r) = z^{-\binom{r}{2}}\sum\limits_{k=0}^n \left(a_k z^{-\binom{k}{2}}\right)z^{\binom{k+r}{2}}}$$

Коефіцієнт біля $x^{n+r}$ у добутку многочленів $A_0(x) = \sum\limits_{k=0}^n a_{n-k}z^{-\binom{n-k}{2}}x^k$ та $A_1(x) = \sum\limits_{k\geq 0}z^{\binom{k}{2}}x^k$ дорівнює $z^{\binom{r}{2}}A(z^r)$. Для обчислення коефіцієнтів $A_0(x)$ та $A_1(x)$ можна використати формулу $z^{\binom{k+1}{2}}=z^{\binom{k}{2}+k}$.

Обчислення в багатьох точках

Припустимо, вам потрібно обчислити $A(x_1), \dots, A(x_n)$. Як було згадано раніше, $A(x) \equiv A(x_i) \pmod{x-x_i}$. Таким чином, можна зробити таке:

1. Побудувати дерево відрізків таке, що у відрізку $[l,r)$ стоїть добуток $P_{l, r}(x) = (x-x_l)(x-x_{l+1})\dots(x-x_{r-1})$.
2. Починаючи з $l=1$ та $r=n+1$ у корені. Нехай $m=\lfloor(l+r)/2\rfloor$. Спускаємося до $[l,m)$ з многочленом $A(x) \pmod{P_{l,m}(x)}$.
3. Це рекурсивно обчислить $A(x_l), \dots, A(x_{m-1})$, тепер зробимо те саме для $[m,r)$ з $A(x) \pmod{P_{m,r}(x)}$.
4. Об'єднати результати першого та другого рекурсивних викликів і повернути їх.

Уся процедура виконається за $O(n \log^2 n)$.

Інтерполяція

Існує пряма формула Лагранжа для інтерполяції многочлена за заданим набором пар $(x_i, y_i)$:

$$\boxed{A(x) = \sum\limits_{i=1}^n y_i \prod\limits_{j \neq i}\dfrac{x-x_j}{x_i - x_j}}$$

Обчислити її напряму важко, але виявляється, що ми можемо обчислити її за $O(n \log^2 n)$ за допомогою підходу «розділяй і володарюй»:

Розглянемо $P(x) = (x-x_1)\dots(x-x_n)$. Щоб знати коефіцієнти знаменників у $A(x)$, ми маємо обчислити добутки на кшталт:

$$P_i = \prod\limits_{j \neq i} (x_i-x_j)$$

Але якщо розглянути похідну $P'(x)$, то виявиться, що $P'(x_i) = P_i$. Таким чином, можна обчислити $P_i$ за допомогою обчислення в точках за $O(n \log^2 n)$.

Тепер розглянемо рекурсивний алгоритм, виконаний на тому самому дереві відрізків, що й в обчисленні в багатьох точках. Він починається в листках зі значенням $\dfrac{y_i}{P_i}$ у кожному листку.

Коли ми повертаємося з рекурсії, ми маємо об'єднати результати з лівої та правої вершин як $A_{l,r} = A_{l,m}P_{m,r} + P_{l,m} A_{m,r}$.

Таким чином, коли ви повернетеся до кореня, ви матимете в ньому саме $A(x)$. Уся процедура також працює за $O(n \log^2 n)$.

НСД та результанти

Припустимо, вам дано многочлени $A(x) = a_0 + a_1 x + \dots + a_n x^n$ та $B(x) = b_0 + b_1 x + \dots + b_m x^m$.

Нехай $\lambda_0, \dots, \lambda_n$ — корені $A(x)$, а $\mu_0, \dots, \mu_m$ — корені $B(x)$, пораховані з урахуванням їхньої кратності.

Ви хочете дізнатися, чи мають $A(x)$ та $B(x)$ якісь спільні корені. Є два взаємопов'язані способи зробити це.

Алгоритм Евкліда

Що ж, у нас уже є стаття про нього. Для довільної області цілісності алгоритм Евкліда можна записати настільки просто:

C++

template<typename T>
T gcd(const T &a, const T &b) {
	return b == T(0) ? a : gcd(b, a % b);
}

Python

# Узагальнений алгоритм Евкліда для будь-якої області цілісності
# (тут T — це многочлен, ціле число тощо; % означає остачу від ділення)
def gcd(a: "T", b: "T") -> "T":
    return a if b == zero_of(b) else gcd(b, a % b)

TypeScript

// Узагальнений алгоритм Евкліда для будь-якої області цілісності;
// % означає остачу від ділення (для многочленів — евклідову остачу)
function gcd<T extends { isZero(): boolean; mod(other: T): T }>(a: T, b: T): T {
  return b.isZero() ? a : gcd(b, a.mod(b));
}

Go

// Узагальнений алгоритм Евкліда для будь-якої області цілісності;
// Mod повертає остачу від ділення (для многочленів — евклідову остачу)
type Euclidean[T any] interface {
	IsZero() bool
	Mod(other T) T
}

func gcd[T Euclidean[T]](a, b T) T {
	if b.IsZero() {
		return a
	}
	return gcd(b, a.Mod(b))
}

Можна довести, що для многочленів $A(x)$ та $B(x)$ він працюватиме за $O(nm)$.

Результант

Обчислимо добуток $A(\mu_0)\cdots A(\mu_m)$. Він дорівнюватиме нулю тоді й лише тоді, коли якесь $\mu_i$ є коренем $A(x)$.

Для симетрії ми також можемо домножити його на $b_m^n$ і переписати весь добуток у такому вигляді:

$$\boxed{\mathcal{R}(A, B) = b_m^n\prod\limits_{j=0}^m A(\mu_j) = b_m^n a_m^n \prod\limits_{i=0}^n \prod\limits_{j=0}^m (\mu_j - \lambda_i)= (-1)^{mn}a_n^m \prod\limits_{i=0}^n B(\lambda_i)}$$

Значення, визначене вище, називається результантом многочленів $A(x)$ та $B(x)$. З означення можна вивести такі властивості:

1. $\mathcal R(A, B) = (-1)^{nm} \mathcal R(B, A)$.
2. $\mathcal R(A, B)= a_n^m b_m^n$, коли $n=0$ або $m=0$.
3. Якщо $b_m=1$, то $\mathcal R(A - CB, B) = \mathcal R(A, B)$ для довільного многочлена $C(x)$ та $n,m \geq 1$.
4. Звідси випливає $\mathcal R(A, B) = b_m^{\deg(A) - \deg(A-CB)}\mathcal R(A - CB, B)$ для довільних $A(x)$, $B(x)$, $C(x)$.

Дивовижним чином це означає, що результант двох многочленів насправді завжди належить тому самому кільцю, що й їхні коефіцієнти!

Також ці властивості дозволяють нам обчислювати результант разом із алгоритмом Евкліда, що працює за $O(nm)$.

C++

template<typename T>
T resultant(poly<T> a, poly<T> b) {
	if(b.is_zero()) {
		return 0;
	} else if(b.deg() == 0) {
		return bpow(b.lead(), a.deg());
	} else {
		int pw = a.deg();
		a %= b;
		pw -= a.deg();
		base mul = bpow(b.lead(), pw) * base((b.deg() & a.deg() & 1) ? -1 : 1);
		base ans = resultant(b, a);
		return ans * mul;
	}
}

Python

# Результант через алгоритм Евкліда; arr — коефіцієнти многочлена
# (arr[i] — коефіцієнт біля x^i), усі операції за модулем MOD
def resultant(a: list[int], b: list[int]) -> int:
    a, b = poly_trim(a), poly_trim(b)
    if poly_deg(b) < 0:  # b порожній (нульовий многочлен)
        return 0
    elif poly_deg(b) == 0:  # b — константа
        return bpow(b[0], poly_deg(a))
    else:
        pw = poly_deg(a)
        a = poly_mod(a, b)  # a %= b
        pw -= poly_deg(a)
        # знак (-1) при непарних степенях обох многочленів
        sign = -1 if (poly_deg(b) & poly_deg(a) & 1) else 1
        mul = bpow(b[-1], pw) * sign % MOD
        return resultant(b, a) * mul % MOD

TypeScript

// Результант через алгоритм Евкліда; arr — коефіцієнти многочлена
// (arr[i] — коефіцієнт біля x^i), усі операції за модулем MOD
function resultant(a: bigint[], b: bigint[]): bigint {
  a = polyTrim(a);
  b = polyTrim(b);
  if (polyDeg(b) < 0) {
    // b порожній (нульовий многочлен)
    return 0n;
  } else if (polyDeg(b) === 0) {
    // b — константа
    return bpow(b[0], polyDeg(a));
  } else {
    let pw = polyDeg(a);
    a = polyMod(a, b); // a %= b
    pw -= polyDeg(a);
    // знак (-1) при непарних степенях обох многочленів
    const sign = (polyDeg(b) & polyDeg(a) & 1) === 1 ? -1n : 1n;
    const mul = ((bpow(b[b.length - 1], pw) * sign) % MOD + MOD) % MOD;
    return ((resultant(b, a) * mul) % MOD + MOD) % MOD;
  }
}

Go

// Результант через алгоритм Евкліда; arr — коефіцієнти многочлена
// (arr[i] — коефіцієнт біля x^i), усі операції за модулем MOD
func resultant(a, b []int64) int64 {
	a, b = polyTrim(a), polyTrim(b)
	if polyDeg(b) < 0 { // b порожній (нульовий многочлен)
		return 0
	} else if polyDeg(b) == 0 { // b — константа
		return bpow(b[0], polyDeg(a))
	}
	pw := polyDeg(a)
	a = polyMod(a, b) // a %= b
	pw -= polyDeg(a)
	// знак (-1) при непарних степенях обох многочленів
	var sign int64 = 1
	if polyDeg(b)&polyDeg(a)&1 == 1 {
		sign = -1
	}
	mul := (bpow(b[len(b)-1], pw)*sign%MOD + MOD) % MOD
	return (resultant(b, a)*mul%MOD + MOD) % MOD
}

Алгоритм Half-GCD

Існує спосіб обчислити НСД та результанти за $O(n \log^2 n)$.

Процедура, що це робить, реалізує лінійне перетворення $2 \times 2$, яке відображає пару многочленів $a(x)$, $b(x)$ в іншу пару $c(x), d(x)$ таку, що $\deg d(x) \leq \frac{\deg a(x)}{2}$. Якщо бути достатньо акуратним, можна обчислити half-GCD будь-якої пари многочленів щонайбільше за $2$ рекурсивні виклики до многочленів, які щонайменше у $2$ рази менші.

Конкретні деталі алгоритму дещо втомливо пояснювати, однак ви можете знайти його реалізацію в бібліотеці як функцію half_gcd.

Після того, як half-GCD реалізовано, ви можете багаторазово застосовувати його до многочленів, доки не зведете їх до пари $\gcd(a, b)$ та $0$.

Задачі

• CodeChef - RNG
• CodeForces - Basis Change
• CodeForces - Permutant
• CodeForces - Medium Hadron Collider

Відеоматеріали