Table of Contents
In this article we will cover common operations that you will probably have to do if you deal with polynomials.
1. Basic Notion and Facts
Consider a polynomial $A(x) = a_0 + a_1 x + \dots + a_n x^n$ such that $a_n \neq 0$.
- For simplicity we will write $A$ instead of $A(x)$ wherever possible, which will be understandable from the context.
- We will define the degree of polynomial $A$ as $\deg A = n$. It is convenient to say that $\deg A = -\infty$ for $A(x) = 0$.
- For arbitrary polynomials $A$ and $B$ it holds that $\deg AB = \deg A + \deg B$.
- Polynomials form an euclidean ring which means that for any polynomials $A$ and $B \neq 0$ we can uniquely represent $A$ as: $$A = D \cdot B + R,~ \deg R < \deg B.$$ Here $R$ is called remainder of $A$ modulo $B$ and $D$ is called quotient.
- If $A$ and $B$ have the same remainder modulo $C$, they’re said to be equivalent modulo $C$, which is denoted as: $$A \equiv B \pmod{C}$$
- For any linear polynomial $x-r$ it holds that: $$A(x) \equiv A(r) \pmod{x-r}$$
- In particular: $$A(r) = 0 \iff A(x) \equiv 0 \pmod {x-r}$$ Which means that $A$ is divisible by $x-r$ $\iff$ $A(r)=0$.
- If $A \equiv B \pmod{C \cdot D}$ then $A \equiv B \pmod{C}$
- $A \equiv a_0 + a_1 x + \dots + a_{k-1} x^{k-1} \pmod{x^k}$
2. Basic implementation
Here you can find the basic implementation of polynomial algebra.
#include <bits/stdc++.h>
using namespace std;
namespace algebra {
const int inf = 1e9;
const int magic = 500; // threshold for sizes to run the naive algo
namespace fft {
const int maxn = 1 << 18;
typedef double ftype;
typedef complex<ftype> point;
point w[maxn];
const ftype pi = acos(-1);
bool initiated = 0;
void init() {
if(!initiated) {
for(int i = 1; i < maxn; i *= 2) {
for(int j = 0; j < i; j++) {
w[i + j] = polar(ftype(1), pi * j / i);
}
}
initiated = 1;
}
}
template<typename T>
void fft(T *in, point *out, int n, int k = 1) {
if(n == 1) {
*out = *in;
} else {
n /= 2;
fft(in, out, n, 2 * k);
fft(in + k, out + n, n, 2 * k);
for(int i = 0; i < n; i++) {
auto t = out[i + n] * w[i + n];
out[i + n] = out[i] - t;
out[i] += t;
}
}
}
template<typename T>
void mul_slow(vector<T> &a, const vector<T> &b) {
vector<T> res(a.size() + b.size() - 1);
for(size_t i = 0; i < a.size(); i++) {
for(size_t j = 0; j < b.size(); j++) {
res[i + j] += a[i] * b[j];
}
}
a = res;
}
template<typename T>
void mul(vector<T> &a, const vector<T> &b) {
if(min(a.size(), b.size()) < magic) {
mul_slow(a, b);
return;
}
init();
static const int shift = 15, mask = (1 << shift) - 1;
size_t n = a.size() + b.size() - 1;
while(__builtin_popcount(n) != 1) {
n++;
}
a.resize(n);
static point A[maxn], B[maxn];
static point C[maxn], D[maxn];
for(size_t i = 0; i < n; i++) {
A[i] = point(a[i] & mask, a[i] >> shift);
if(i < b.size()) {
B[i] = point(b[i] & mask, b[i] >> shift);
} else {
B[i] = 0;
}
}
fft(A, C, n); fft(B, D, n);
for(size_t i = 0; i < n; i++) {
point c0 = C[i] + conj(C[(n - i) % n]);
point c1 = C[i] - conj(C[(n - i) % n]);
point d0 = D[i] + conj(D[(n - i) % n]);
point d1 = D[i] - conj(D[(n - i) % n]);
A[i] = c0 * d0 - point(0, 1) * c1 * d1;
B[i] = c0 * d1 + d0 * c1;
}
fft(A, C, n); fft(B, D, n);
reverse(C + 1, C + n);
reverse(D + 1, D + n);
int t = 4 * n;
for(size_t i = 0; i < n; i++) {
int64_t A0 = llround(real(C[i]) / t);
T A1 = llround(imag(D[i]) / t);
T A2 = llround(imag(C[i]) / t);
a[i] = A0 + (A1 << shift) + (A2 << 2 * shift);
}
return;
}
}
template<typename T>
T bpow(T x, size_t n) {
return n ? n % 2 ? x * bpow(x, n - 1) : bpow(x * x, n / 2) : T(1);
}
template<typename T>
T bpow(T x, size_t n, T m) {
return n ? n % 2 ? x * bpow(x, n - 1, m) % m : bpow(x * x % m, n / 2, m) : T(1);
}
template<typename T>
T gcd(const T &a, const T &b) {
return b == T(0) ? a : gcd(b, a % b);
}
template<typename T>
T nCr(T n, int r) { // runs in O(r)
T res(1);
for(int i = 0; i < r; i++) {
res *= (n - T(i));
res /= (i + 1);
}
return res;
}
template<int m>
struct modular {
int64_t r;
modular() : r(0) {}
modular(int64_t rr) : r(rr) {if(abs(r) >= m) r %= m; if(r < 0) r += m;}
modular inv() const {return bpow(*this, m - 2);}
modular operator * (const modular &t) const {return (r * t.r) % m;}
modular operator / (const modular &t) const {return *this * t.inv();}
modular operator += (const modular &t) {r += t.r; if(r >= m) r -= m; return *this;}
modular operator -= (const modular &t) {r -= t.r; if(r < 0) r += m; return *this;}
modular operator + (const modular &t) const {return modular(*this) += t;}
modular operator - (const modular &t) const {return modular(*this) -= t;}
modular operator *= (const modular &t) {return *this = *this * t;}
modular operator /= (const modular &t) {return *this = *this / t;}
bool operator == (const modular &t) const {return r == t.r;}
bool operator != (const modular &t) const {return r != t.r;}
operator int64_t() const {return r;}
};
template<int T>
istream& operator >> (istream &in, modular<T> &x) {
return in >> x.r;
}
template<typename T>
struct poly {
vector<T> a;
void normalize() { // get rid of leading zeroes
while(!a.empty() && a.back() == T(0)) {
a.pop_back();
}
}
poly(){}
poly(T a0) : a{a0}{normalize();}
poly(vector<T> t) : a(t){normalize();}
poly operator += (const poly &t) {
a.resize(max(a.size(), t.a.size()));
for(size_t i = 0; i < t.a.size(); i++) {
a[i] += t.a[i];
}
normalize();
return *this;
}
poly operator -= (const poly &t) {
a.resize(max(a.size(), t.a.size()));
for(size_t i = 0; i < t.a.size(); i++) {
a[i] -= t.a[i];
}
normalize();
return *this;
}
poly operator + (const poly &t) const {return poly(*this) += t;}
poly operator - (const poly &t) const {return poly(*this) -= t;}
poly mod_xk(size_t k) const { // get same polynomial mod x^k
k = min(k, a.size());
return vector<T>(begin(a), begin(a) + k);
}
poly mul_xk(size_t k) const { // multiply by x^k
poly res(*this);
res.a.insert(begin(res.a), k, 0);
return res;
}
poly div_xk(size_t k) const { // divide by x^k, dropping coefficients
k = min(k, a.size());
return vector<T>(begin(a) + k, end(a));
}
poly substr(size_t l, size_t r) const { // return mod_xk(r).div_xk(l)
l = min(l, a.size());
r = min(r, a.size());
return vector<T>(begin(a) + l, begin(a) + r);
}
poly inv(size_t n) const { // get inverse series mod x^n
assert(!is_zero());
poly ans = a[0].inv();
size_t a = 1;
while(a < n) {
poly C = (ans * mod_xk(2 * a)).substr(a, 2 * a);
ans -= (ans * C).mod_xk(a).mul_xk(a);
a *= 2;
}
return ans.mod_xk(n);
}
poly operator *= (const poly &t) {fft::mul(a, t.a); normalize(); return *this;}
poly operator * (const poly &t) const {return poly(*this) *= t;}
poly reverse(size_t n, bool rev = 0) const { // reverses and leaves only n terms
poly res(*this);
if(rev) { // If rev = 1 then tail goes to head
res.a.resize(max(n, res.a.size()));
}
std::reverse(res.a.begin(), res.a.end());
return res.mod_xk(n);
}
pair<poly, poly> divmod_slow(const poly &b) const { // when divisor or quotient is small
vector<T> A(a);
vector<T> res;
while(A.size() >= b.a.size()) {
res.push_back(A.back() / b.a.back());
if(res.back() != T(0)) {
for(size_t i = 0; i < b.a.size(); i++) {
A[A.size() - i - 1] -= res.back() * b.a[b.a.size() - i - 1];
}
}
A.pop_back();
}
std::reverse(begin(res), end(res));
return {res, A};
}
pair<poly, poly> divmod(const poly &b) const { // returns quotiend and remainder of a mod b
if(deg() < b.deg()) {
return {poly{0}, *this};
}
int d = deg() - b.deg();
if(min(d, b.deg()) < magic) {
return divmod_slow(b);
}
poly D = (reverse(d + 1) * b.reverse(d + 1).inv(d + 1)).mod_xk(d + 1).reverse(d + 1, 1);
return {D, *this - D * b};
}
poly operator / (const poly &t) const {return divmod(t).first;}
poly operator % (const poly &t) const {return divmod(t).second;}
poly operator /= (const poly &t) {return *this = divmod(t).first;}
poly operator %= (const poly &t) {return *this = divmod(t).second;}
poly operator *= (const T &x) {
for(auto &it: a) {
it *= x;
}
normalize();
return *this;
}
poly operator /= (const T &x) {
for(auto &it: a) {
it /= x;
}
normalize();
return *this;
}
poly operator * (const T &x) const {return poly(*this) *= x;}
poly operator / (const T &x) const {return poly(*this) /= x;}
void print() const {
for(auto it: a) {
cout << it << ' ';
}
cout << endl;
}
T eval(T x) const { // evaluates in single point x
T res(0);
for(int i = int(a.size()) - 1; i >= 0; i--) {
res *= x;
res += a[i];
}
return res;
}
T& lead() { // leading coefficient
return a.back();
}
int deg() const { // degree
return a.empty() ? -inf : a.size() - 1;
}
bool is_zero() const { // is polynomial zero
return a.empty();
}
T operator [](int idx) const {
return idx >= (int)a.size() || idx < 0 ? T(0) : a[idx];
}
T& coef(size_t idx) { // mutable reference at coefficient
return a[idx];
}
bool operator == (const poly &t) const {return a == t.a;}
bool operator != (const poly &t) const {return a != t.a;}
poly deriv() { // calculate derivative
vector<T> res;
for(int i = 1; i <= deg(); i++) {
res.push_back(T(i) * a[i]);
}
return res;
}
poly integr() { // calculate integral with C = 0
vector<T> res = {0};
for(int i = 0; i <= deg(); i++) {
res.push_back(a[i] / T(i + 1));
}
return res;
}
size_t leading_xk() const { // Let p(x) = x^k * t(x), return k
if(is_zero()) {
return inf;
}
int res = 0;
while(a[res] == T(0)) {
res++;
}
return res;
}
poly log(size_t n) { // calculate log p(x) mod x^n
assert(a[0] == T(1));
return (deriv().mod_xk(n) * inv(n)).integr().mod_xk(n);
}
poly exp(size_t n) { // calculate exp p(x) mod x^n
if(is_zero()) {
return T(1);
}
assert(a[0] == T(0));
poly ans = T(1);
size_t a = 1;
while(a < n) {
poly C = ans.log(2 * a).div_xk(a) - substr(a, 2 * a);
ans -= (ans * C).mod_xk(a).mul_xk(a);
a *= 2;
}
return ans.mod_xk(n);
}
poly pow_slow(size_t k, size_t n) { // if k is small
return k ? k % 2 ? (*this * pow_slow(k - 1, n)).mod_xk(n) : (*this * *this).mod_xk(n).pow_slow(k / 2, n) : T(1);
}
poly pow(size_t k, size_t n) { // calculate p^k(n) mod x^n
if(is_zero()) {
return *this;
}
if(k < magic) {
return pow_slow(k, n);
}
int i = leading_xk();
T j = a[i];
poly t = div_xk(i) / j;
return bpow(j, k) * (t.log(n) * T(k)).exp(n).mul_xk(i * k).mod_xk(n);
}
poly mulx(T x) { // component-wise multiplication with x^k
T cur = 1;
poly res(*this);
for(int i = 0; i <= deg(); i++) {
res.coef(i) *= cur;
cur *= x;
}
return res;
}
poly mulx_sq(T x) { // component-wise multiplication with x^{k^2}
T cur = x;
T total = 1;
T xx = x * x;
poly res(*this);
for(int i = 0; i <= deg(); i++) {
res.coef(i) *= total;
total *= cur;
cur *= xx;
}
return res;
}
vector<T> chirpz_even(T z, int n) { // P(1), P(z^2), P(z^4), ..., P(z^2(n-1))
int m = deg();
if(is_zero()) {
return vector<T>(n, 0);
}
vector<T> vv(m + n);
T zi = z.inv();
T zz = zi * zi;
T cur = zi;
T total = 1;
for(int i = 0; i <= max(n - 1, m); i++) {
if(i <= m) {vv[m - i] = total;}
if(i < n) {vv[m + i] = total;}
total *= cur;
cur *= zz;
}
poly w = (mulx_sq(z) * vv).substr(m, m + n).mulx_sq(z);
vector<T> res(n);
for(int i = 0; i < n; i++) {
res[i] = w[i];
}
return res;
}
vector<T> chirpz(T z, int n) { // P(1), P(z), P(z^2), ..., P(z^(n-1))
auto even = chirpz_even(z, (n + 1) / 2);
auto odd = mulx(z).chirpz_even(z, n / 2);
vector<T> ans(n);
for(int i = 0; i < n / 2; i++) {
ans[2 * i] = even[i];
ans[2 * i + 1] = odd[i];
}
if(n % 2 == 1) {
ans[n - 1] = even.back();
}
return ans;
}
template<typename iter>
vector<T> eval(vector<poly> &tree, int v, iter l, iter r) { // auxiliary evaluation function
if(r - l == 1) {
return {eval(*l)};
} else {
auto m = l + (r - l) / 2;
auto A = (*this % tree[2 * v]).eval(tree, 2 * v, l, m);
auto B = (*this % tree[2 * v + 1]).eval(tree, 2 * v + 1, m, r);
A.insert(end(A), begin(B), end(B));
return A;
}
}
vector<T> eval(vector<T> x) { // evaluate polynomial in (x1, ..., xn)
int n = x.size();
if(is_zero()) {
return vector<T>(n, T(0));
}
vector<poly> tree(4 * n);
build(tree, 1, begin(x), end(x));
return eval(tree, 1, begin(x), end(x));
}
template<typename iter>
poly inter(vector<poly> &tree, int v, iter l, iter r, iter ly, iter ry) { // auxiliary interpolation function
if(r - l == 1) {
return {*ly / a[0]};
} else {
auto m = l + (r - l) / 2;
auto my = ly + (ry - ly) / 2;
auto A = (*this % tree[2 * v]).inter(tree, 2 * v, l, m, ly, my);
auto B = (*this % tree[2 * v + 1]).inter(tree, 2 * v + 1, m, r, my, ry);
return A * tree[2 * v + 1] + B * tree[2 * v];
}
}
};
template<typename T>
poly<T> operator * (const T& a, const poly<T>& b) {
return b * a;
}
template<typename T>
poly<T> xk(int k) { // return x^k
return poly<T>{1}.mul_xk(k);
}
template<typename T>
T resultant(poly<T> a, poly<T> b) { // computes resultant of a and 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();
T mul = bpow(b.lead(), pw) * T((b.deg() & a.deg() & 1) ? -1 : 1);
T ans = resultant(b, a);
return ans * mul;
}
}
template<typename iter>
poly<typename iter::value_type> kmul(iter L, iter R) { // computes (x-a1)(x-a2)...(x-an) without building tree
if(R - L == 1) {
return vector<typename iter::value_type>{-*L, 1};
} else {
iter M = L + (R - L) / 2;
return kmul(L, M) * kmul(M, R);
}
}
template<typename T, typename iter>
poly<T> build(vector<poly<T>> &res, int v, iter L, iter R) { // builds evaluation tree for (x-a1)(x-a2)...(x-an)
if(R - L == 1) {
return res[v] = vector<T>{-*L, 1};
} else {
iter M = L + (R - L) / 2;
return res[v] = build(res, 2 * v, L, M) * build(res, 2 * v + 1, M, R);
}
}
template<typename T>
poly<T> inter(vector<T> x, vector<T> y) { // interpolates minimum polynomial from (xi, yi) pairs
int n = x.size();
vector<poly<T>> tree(4 * n);
return build(tree, 1, begin(x), end(x)).deriv().inter(tree, 1, begin(x), end(x), begin(y), end(y));
}
};
using namespace algebra;
const int mod = 1e9 + 7;
typedef modular<mod> base;
typedef poly<base> polyn;
using namespace algebra;
signed main() {
ios::sync_with_stdio(0);
cin.tie(0);
int n = 100000;
polyn a;
vector<base> x;
for(int i = 0; i <= n; i++) {
a.a.push_back(1 + rand() % 100);
x.push_back(1 + rand() % (2 * n));
}
sort(begin(x), end(x));
x.erase(unique(begin(x), end(x)), end(x));
auto b = a.eval(x);
cout << clock() / double(CLOCKS_PER_SEC) << endl;
auto c = inter(x, b);
polyn md = kmul(begin(x), end(x));
cout << clock() / double(CLOCKS_PER_SEC) << endl;
assert(c == a % md);
return 0;
}
It supports all trivial operations and some other useful methods. The main class is poly for polynomials with coefficients of class T.
All arithmetic operation +, -, *, % and / are supported, % and / standing for remainder and quotient in integer division.
There is also the class modular for performing arithmetic operations on remainders modulo a prime number m.
Other useful functions:
- deriv(): computes the derivative $P'(x)$ of $P(x)$.
- integr(): computes the indefinite integral $Q(x) = \int P(x)$ of $P(x)$ such that $Q(0)=0$.
- inv(size_t n): calculate the first $n$ coefficients of $P^{-1}(x)$ in $O(n \log n)$.
- log(size_t n): calculate the first $n$ coefficients of $\ln P(x)$ in $O(n \log n)$.
- exp(size_t n): calculate the first $n$ coefficients of $\exp P(x)$ in $O(n \log n)$.
- pow(size_t k, size_t n): calculate the first $n$ coefficients for $P^{k}(x)$ in $O(n \log nk)$.
- deg(): returns the degree of $P(x)$.
- lead(): returns the coefficient of $x^{\deg P(x)}$.
- resultant(poly a, poly b): computes the resultant of $a$ and $b$ in $O(|a| \cdot |b|)$.
- bpow(T x, size_t n): computes $x^n$.
- bpow(T x, size_t n, T m): computes $x^n \pmod{m}$.
- chirpz(T z, size_t n): computes $P(1), P(z), P(z^2), \dots, P(z^{n-1})$ in $O(n \log n)$.
- vector eval(vector x): evaluates $P(x_1), \dots, P(x_n)$ in $O(n \log^2 n)$.
- poly inter(vector x, vector y): interpolates a polynomial by a set of pairs $P(x_i) = y_i$ in $O(n \log^2 n)$.
- And some more, feel free to explore the code!
3. Arithmetic
3.1. Multiplication
The very core operation is the multiplication of two polynomials, that is, given polynomial $A$ and $B$: $$A = a_0 + a_1 x + \dots + a_n x^n$$ $$B = b_0 + b_1 x + \dots + b_m x^m$$ You have to compute polynomial $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}$$ It can be computed in $O(n \log n)$ via the Fast Fourier transform and almost all methods here will use it as subroutine.
3.2. Inverse series
If $A(0) \neq 0$ there always exists an infinite series $A^{-1}(x) = \sum\limits_{i=0}^\infty a_i’x^i$ such that $A^{-1} A = 1$.
It may be reasonable for us to calculate first $k$ coefficients of $A^{-1}$:
- Let’s say that $A^{-1} \equiv B_k \pmod{x^{a}}$. That means that $A B_k \equiv 1 \pmod {x^{a}}$.
- We want to find $B_{k+1} \equiv B_k + x^{a}C \pmod{x^{2a}}$ such that $A B_{k+1} \equiv 1 \pmod{x^{2a}}$: $$A(B_k + x^{a}C) \equiv 1 \pmod{x^{2a}}$$
- Note that since $A B_k \equiv 1 \pmod{x^{a}}$ it also holds that $A B_k \equiv 1 + x^a D \pmod{x^{2a}}$. Thus: $$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}$$
- From this we obtain that: $$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}}}$$
Thus starting with $B_0 \equiv a_0^{-1} \pmod x$ we will compute the sequence $B_k$ such that $AB_k \equiv 1 \pmod{x^{2^k}}$ with the complexity: $$T(n) = T(n/2) + O(n \log n) = O(n \log n)$$
3.3. Division with remainder
Consider two polynomials $A(x)$ and $B(x)$ of degrees $n$ and $m$. As it was said earlier you can rewrite $A(x)$ as:
$$A(x) = B(x) D(x) + R(x), \deg R < \deg B$$
Let $n \geq m$, then you can immediately find out that $\deg D = n – m$ and that leading $n-m+1$ coefficients of $A$ don’t influence $R$.
That means that you can recover $D(x)$ from the largest $n-m+1$ coefficients of $A(x)$ and $B(x)$ if you consider it as a system of equations.
The formal way to do it is to consider the reversed polynomials: $$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}$$
Using these terms you can rewrite that statement about the largest $n-m+1$ coefficients as:
$$A^R(x) \equiv B^R(x) D^R(x) \pmod{x^{n-m+1}}$$
From which you can unambiguously recover all coefficients of $D(x)$:
$$\boxed{D^R(x) \equiv A^R(x) (B^R(x))^{-1} \pmod{x^{n-m+1}}}$$
And from this in turn you can easily recover $R(x)$ as $R(x) = A(x) – B(x)D(x)$.
4. Calculating functions of polynomial
4.1. Newton’s method
Let’s generalize the inverse series approach. You want to find a polynomial $P(x)$ satisfying $F(P) = 0$ where $F(x)$ is some function represented as: $$F(x) = \sum\limits_{i=0}^\infty \alpha_i (x-\beta)^k$$
Where $\beta$ is some constant. It can be proven that if we introduce a new formal variable $y$, we can express $F(x)$ as: $$F(x) = F(y) + (x-y)F'(y) + (x-y)^2 G(x,y)$$
Where $F'(x)$ is the derivative formal power series defined as $F'(x) = \sum\limits_{i=0}^\infty (k+1)\alpha_{i+1}(x-\beta)^k$ and $G(x, y)$ is some formal power series of $x$ and $y$.
Given this we can find the coefficients of the solution iteratively:
- Assume that $F(Q_k) \equiv 0 \pmod{x^{a}}$, we want to find $Q_{k+1} \equiv Q_k + x^a C \pmod{x^{2a}}$ such that $F(Q_{k+1}) \equiv 0 \pmod{x^{2a}}$.
- Pasting $x = Q_{k+1}$ and $y=Q_k$ in the formula above we get: $$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}$$
- Since $Q_{k+1} – Q_k \equiv 0 \pmod{x^a}$ we can say that $(Q_{k+1} – Q_k)^2 \equiv 0 \pmod{x^{2a}}$, thus: $$0 \equiv F(Q_{k+1}) \equiv F(Q_k) + (Q_{k+1} – Q_k) F'(Q_k) \pmod{x^{2a}}$$
- From the last formula we derive the value of $Q_{k+1}$: $$\boxed{Q_{k+1} = Q_k – \dfrac{F(Q_k)}{F'(Q_k)} \pmod{x^{2a}}}$$
Thus knowing how to invert arbitrary polynomial and how to compute $F(Q_k)$ quickly, we can find $n$ coefficients of $P$ with the complexity: $$T(n) = T(n/2) + f(n)$$ Where $f(n)$ is the maximum of $O(n \log n)$ needed to invert series and time needed to compute $F(Q_k)$ which is usually also $O(n \log n)$.
4.2. Logarithm
For the function $\ln P(x)$ it’s known that: $$\boxed{(\ln P(x))’ = \dfrac{P'(x)}{P(x)}}$$ Thus we can calculate $n$ coefficients of $\ln P(x)$ in $O(n \log n)$.
4.3. Inverse series
Turns out, we can get the formula for $A^{-1}$ using Newton’s method. For this we take the equation $A=Q^{-1}$, thus: $$F(Q) = Q^{-1} – A$$ $$F'(Q) = -Q^{-2}$$ $$\boxed{Q_{k+1} \equiv Q_k(2-AQ_k) \pmod{x^{2^{k+1}}}}$$
4.4. Exponent
Let’s learn to calculate $e^{P(x)}=Q(x)$. It should hold that $\ln Q = P$, thus: $$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}}}}$$
4.5. $k$-th power
Now we need to calculate $P^k(x)=Q$. This may be done via the following formula: $$Q = \exp\left[k \ln P(x)\right]$$ Note though, that you can calculate the logarithms and the exponents correctly only if you can find some initial $Q_0$.
To find it, you should calculate the logarithm or the exponent of the constant coefficient of the polynomial.
But the only reasonable way to do it is if $P(0)=1$ for $Q = \ln P$ so $Q(0)=0$ and if $P(0)=0$ for $Q = e^P$ so $Q(0)=1$.
Thus you can use formula above only if $P(0) = 1$. Otherwise if $P(x) = \alpha x^t T(x)$ where $T(0)=1$ you can write that: $$\boxed{P^k(x) = \alpha^kx^{kt} \exp[k \ln T(x)]}$$ Note that you also can calculate some $k$-th root of a polynomial if you can calculate $\sqrt[k]{\alpha}$, for example for $\alpha=1$.
5. Evaluation and Interpolation
5.1. Chirp-z Transform
For the particular case when you need to evaluate a polynomial in the points $x_r = z^{2r}$ you can do the following:
$$A(z^{2r}) = \sum\limits_{k=0}^n a_k z^{2kr}$$
Let’s substitute $2kr = r^2+k^2-(r-k)^2$. Then this sum rewrites as:
$$\boxed{A(z^{2r}) = z^{r^2}\sum\limits_{k=0}^n (a_k z^{k^2}) z^{-(r-k)^2}}$$
Which is up to the factor $z^{r^2}$ equal to the convolution of the sequences $u_k = a_k z^{k^2}$ and $v_k = z^{-k^2}$.
Note that $u_k$ has indexes from $0$ to $n$ here and $v_k$ has indexes from $-n$ to $m$ where $m$ is the maximum power of $z$ which you need.
Now if you need to evaluate a polynomial in the points $x_r = z^{2r+1}$ you can reduce it to the previous task by the transformation $a_k \to a_k z^k$.
It gives us an $O(n \log n)$ algorithm when you need to compute values in powers of $z$, thus you may compute the DFT for non-powers of two.
5.2. Multi-point Evaluation
Assume you need to calculate $A(x_1), \dots, A(x_n)$. As mentioned earlier, $A(x) \equiv A(x_i) \pmod{x-x_i}$. Thus you may do the following:
- Compute a segment tree such that in the segment $[l,r)$ stands the product $P_{l, r}(x) = (x-x_l)(x-x_{l+1})\dots(x-x_{r-1})$.
- Starting with $l=1$ and $r=n$ at the root node. Let $m=\lfloor(l+r)/2\rfloor$. Let’s move down to $[l,m)$ with the polynomial $A(x) \pmod{P_{l,m}(x)}$.
- This will recursively compute $A(x_l), \dots, A(x_{m-1})$, now do the same for $[m,r)$ with $A(x) \pmod{P_{m,r}(x)}$.
- Concatenate the results from the first and second recursive call and return them.
The whole procedure will run in $O(n \log^2 n)$.
5.3. Interpolation
There’s a direct formula by Lagrange to interpolate a polynomial, given set of pairs $(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}}$$
Computing it directly is a hard thing but turns out, we may compute it in $O(n \log^2 n)$ with a divide and conquer approach:
Consider $P(x) = (x-x_1)\dots(x-x_n)$. To know the coefficients of the denominators in $A(x)$ we should compute products like: $$P_i = \prod\limits_{j \neq i} (x_i-x_j)$$
But if you consider the derivative $P'(x)$ you’ll find out that $P'(x_i) = P_i$. Thus you can compute $P_i$’s via evaluation in $O(n \log^2 n)$.
Now consider the recursive algorithm done on same segment tree as in the multi-point evaluation. It starts in the leaves with the value $\dfrac{y_i}{P_i}$ in each leaf.
When we return from the recursion we should merge the results from the left and the right vertices as $A_{l,r} = A_{l,m}P_{m,r} + P_{l,m} A_{m,r}$.
In this way when you return back to the root you’ll have exactly $A(x)$ in it. The total procedure also works in $O(n \log^2 n)$.
6. GCD and Resultants
Assume you’re given polynomials $A(x) = a_0 + a_1 x + \dots + a_n x^n$ and $B(x) = b_0 + b_1 x + \dots + b_m x^m$.
Let $\lambda_0, \dots, \lambda_n$ be the roots of $A(x)$ and let $\mu_0, \dots, \mu_m$ be the roots of $B(x)$ counted with their multiplicities.
You want to know if $A(x)$ and $B(x)$ have any roots in common. There are two interconnected ways to do that.
6.1. Euclidean algorithm
Well, we already have an article about it. For an arbitrary euclidean domain you can write the Euclidean algorithm as easy as:
template<typename T>
T gcd(const T &a, const T &b) {
return b == T(0) ? a : gcd(b, a % b);
}
It can be proven that for polynomials $A(x)$ and $B(x)$ it will work in $O(nm)$.
6.2. Resultant
Let’s calculate the product $A(\mu_0)\cdots A(\mu_m)$. It will be equal to zero if and only if some $\mu_i$ is the root of $A(x)$.
For symmetry we can also multiply it with $b_m^n$ and rewrite the whole product in the following form: $$\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)}$$
The value defined above is called the resultant of the polynomials $A(x)$ and $B(x)$. From the definition you may find the following properties:
- $\mathcal R(A, B) = (-1)^{nm} \mathcal R(B, A)$.
- $\mathcal R(A, B)= a_n^m b_m^n$ when $n=0$ or $m=0$.
- If $b_m=1$ then $\mathcal R(A – CB, B) = \mathcal R(A, B)$ for an arbitrary polynomial $C(x)$ and $n,m \geq 1$.
- From this follows $\mathcal R(A, B) = b_m^{\deg(A) – \deg(A-CB)}\mathcal R(A – CB, B)$ for arbitrary $A(x)$, $B(x)$, $C(x)$.
Miraculously it means that resultant of two polynomials is actually always from the same ring as their coefficients!
Also these properties allow us to calculate the resultant alongside the Euclidean algorithm, which works in $O(nm)$.
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;
}
}
6.3. Half-GCD algorithm
There is a way to calculate the GCD and resultants in $O(n \log^2 n)$. To do this you should note that if you consider $a(x) = a_0 + x^k a_1$ and $b(x) = b_0 + x^k b_1$ where $k=\min(\deg a, \deg b)/2$ then basically the first few operations of Euclidean algorithm on $a(x)$ and $b(x)$ are defined by the Euclidean algorithm on $a_1(x)$ and $b_1(x)$ for which you may also calculate GCD recursively and then somehow memorize linear transforms you made with them and apply it to $a(x)$ and $b(x)$ to lower the degrees of polynomials. Implementation of this algorithm seems pretty tedious and technical thus it’s not considered in this article yet.
