This page covers common elementary operations on polynomials. Specifically, it includes the following:
- Polynomial inversion
- Polynomial square roots
- Polynomial division
- Polynomial remainder
- Polynomial exponentials
- Polynomial logarithms
- Polynomial trigonometric functions
- Polynomial inverse trigonometric functions
Elementary and non-elementary functions
The definition of an elementary function is as follows1:
If there exists a mapping in a field satisfying:
then this field is called a differential field.
If a function on the differential field satisfies any one of the following conditions, then is called an elementary function:
- is an algebraic function over .
- is an exponential function over , meaning that there exists such that .
- is a logarithmic function over , meaning that there exists such that .
The following are common elementary functions:
-
Algebraic functions: functions for which there exists a finite-degree polynomial such that , such as ,,,.
-
Exponential functions
-
Logarithmic functions
-
Trigonometric functions
-
Inverse trigonometric functions
-
Hyperbolic functions
-
Inverse hyperbolic functions
-
Compositions of the above functions, such as:
The following is a common non-elementary function:
-
The error function:
Polynomial Inversion
Given a polynomial , find .
Solution
Doubling Method
First, it is easy to see that
Assume that the inverse of modulo has already been found. Then:
Squaring both sides gives:
Multiplying both sides by and rearranging gives:
The result follows by recursion.
Time complexity
Newton's Method
See Newton's Method.
Graeffe's Method
To find , consider
It is enough to find to recover , because is an even function. The time complexity is the same as above.
Code
Polynomial Inversion
constexpr int MAXN = 262144;
constexpr int mod = 998244353;
using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;
void polyinv(const poly &h, const int n, poly &f) {
/* f = 1 / h = f_0 (2 - f_0 h) */
static poly_t inv_t;
std::fill(f, f + n + n, 0);
f[0] = fpow(h[0], mod - 2);
for (int t = 2; t <= n; t <<= 1) {
const int t2 = t << 1;
std::copy(h, h + t, inv_t);
std::fill(inv_t + t, inv_t + t2, 0);
DFT(f, t2);
DFT(inv_t, t2);
for (int i = 0; i != t2; ++i)
f[i] = (i64)f[i] * sub(2, (i64)f[i] * inv_t[i] % mod) % mod;
IDFT(f, t2);
std::fill(f + t, f + t2, 0);
}
}
Example
- Counting labeled simple undirected connected graphs: 「POJ 1737」Connected Graph
Polynomial Square Roots
Given a polynomial , find satisfying:
Solution
Doubling Method
First consider the case where is not .
It is easy to see that:
If has no square root, then the polynomial has no square root.
may have multiple square roots; choosing different roots yields different .
Assume that the square root of modulo has already been found. Then:
Compute by doubling.
Time complexity
Another version with a smaller constant is to maintain alongside during doubling, instead of computing the inverse each time.
When , a quadratic-residue algorithm may be needed to compute .
The method above requires knowing the inverse of , so its constant term cannot be .
If , factor as , where .
-
If is odd, then has no square root.
-
If is even, find the square root of , and then obtain .
Luogu Template Problem P5205 【Template】Polynomial Square Root Reference Code
#include <algorithm>
#include <iostream>
using namespace std;
constexpr int MAXN = 1 << 20, mod = 998244353;
int a[MAXN], b[MAXN], g[MAXN], gg[MAXN];
int qpow(int x, int y) { // fast exponentiation
int ans = 1;
while (y) {
if (y & 1) {
ans = (long long)1 * ans * x % mod;
}
x = (long long)1 * x * x % mod;
y >>= 1;
}
return ans;
}
int inv2 = qpow(2, mod - 2); // modular inverse
void change(int *f, int len) {
for (int i = 1, j = len / 2; i < len - 1; i++) {
if (i < j) {
swap(f[i], f[j]);
}
int k = len / 2;
while (j >= k) {
j -= k;
k /= 2;
}
if (j < k) {
j += k;
}
}
}
void NTT(int *f, int len, int type) { // NTT
change(f, len);
for (int q = 2; q <= len; q <<= 1) {
int nxt = qpow(3, (mod - 1) / q);
for (int i = 0; i < len; i += q) {
int w = 1;
for (int k = i; k < i + (q >> 1); k++) {
int x = f[k];
int y = (long long)1 * w * f[k + (q / 2)] % mod;
f[k] = (x + y) % mod;
f[k + (q / 2)] = (x - y + mod) % mod;
w = (long long)1 * w * nxt % mod;
}
}
}
if (type == -1) {
reverse(f + 1, f + len);
int iv = qpow(len, mod - 2);
for (int i = 0; i < len; i++) {
f[i] = (long long)1 * f[i] * iv % mod;
}
}
}
void inv(int deg, int *f, int *h) { // polynomial inversion
if (deg == 1) {
h[0] = qpow(f[0], mod - 2);
return;
}
inv((deg + 1) >> 1, f, h);
int len = 1;
while (len < deg * 2) { // doubling
len *= 2;
}
copy(f, f + deg, gg);
fill(gg + deg, gg + len, 0);
NTT(gg, len, 1);
NTT(h, len, 1);
for (int i = 0; i < len; i++) {
h[i] = (long long)1 * (2 - (long long)1 * gg[i] * h[i] % mod + mod) % mod *
h[i] % mod;
}
NTT(h, len, -1);
fill(h + deg, h + len, 0);
}
int n, t[MAXN];
// deg: degree
// f: input array whose square root is taken
// h: output array
void sqrt(int deg, int *f, int *h) {
if (deg == 1) {
h[0] = 1;
return;
}
sqrt((deg + 1) >> 1, f, h);
int len = 1;
while (len < deg * 2) { // doubling
len *= 2;
}
fill(g, g + len, 0);
inv(deg, h, g);
copy(f, f + deg, t);
fill(t + deg, t + len, 0);
NTT(t, len, 1);
NTT(g, len, 1);
NTT(h, len, 1);
for (int i = 0; i < len; i++) {
h[i] = (long long)1 * inv2 *
((long long)1 * h[i] % mod + (long long)1 * g[i] * t[i] % mod) % mod;
}
NTT(h, len, -1);
fill(h + deg, h + len, 0);
}
int main() {
cin >> n;
for (int i = 0; i < n; i++) {
cin >> a[i];
}
sqrt(n, a, b);
for (int i = 0; i < n; i++) {
cout << b[i] << ' ';
}
return 0;
}
Newton's Method
See Newton's Method.
Example
Polynomial Division & Modulo
Given polynomials , find the quotient and remainder when is divided by .
Solution
If the effect of can be eliminated, the problem can be solved directly using polynomial inversion.
Consider constructing the transformation
It can be seen that this simply reverses the coefficients of .
Let .
Substitute for in and multiply both sides by to obtain:
Notice that the coefficients of in the equation above are multiplied by , so taking the equation modulo eliminates the effect of .
Also, since the degree of is , is unaffected.
Therefore:
The quotient can be obtained using polynomial inversion, and substituting it back gives .
Time complexity .
Polynomial Logarithms & Exponentials
Given a polynomial , find and modulo .
Solution
Ordinary Method
Polynomial logarithm
First, for a polynomial , if exists, then by its definition, it must satisfy:
Differentiate and then integrate to obtain:
Polynomial differentiation and integration each take time, while inversion takes time. Thus, computing the polynomial takes time.
Polynomial exponential
First, for a polynomial , if exists, then it must satisfy:
Otherwise, the constant term of does not converge.
Differentiate to obtain:
Comparing coefficients on both sides gives:
This can be solved using divide-and-conquer FFT.
Time complexity .
Newton's Method
Using Newton's Method, the polynomial can be computed in time.
Code
Polynomial ln/exp
constexpr int MAXN = 262144;
constexpr int mod = 998244353;
using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;
void derivative(const poly &h, const int n, poly &f) {
for (int i = 1; i != n; ++i) f[i - 1] = (i64)h[i] * i % mod;
f[n - 1] = 0;
}
void integrate(const poly &h, const int n, poly &f) {
for (int i = n - 1; i; --i) f[i] = (i64)h[i - 1] * inv[i] % mod;
f[0] = 0; /* C */
}
void polyln(const poly &h, const int n, poly &f) {
/* f = ln h = ∫ h' / h dx */
assert(h[0] == 1);
static poly_t ln_t;
const int t = n << 1;
derivative(h, n, ln_t);
std::fill(ln_t + n, ln_t + t, 0);
polyinv(h, n, f);
DFT(ln_t, t);
DFT(f, t);
for (int i = 0; i != t; ++i) ln_t[i] = (i64)ln_t[i] * f[i] % mod;
IDFT(ln_t, t);
integrate(ln_t, n, f);
}
void polyexp(const poly &h, const int n, poly &f) {
/* f = exp(h) = f_0 (1 - ln f_0 + h) */
assert(h[0] == 0);
static poly_t exp_t;
std::fill(f, f + n + n, 0);
f[0] = 1;
for (int t = 2; t <= n; t <<= 1) {
const int t2 = t << 1;
polyln(f, t, exp_t);
exp_t[0] = sub(pls(h[0], 1), exp_t[0]);
for (int i = 1; i != t; ++i) exp_t[i] = sub(h[i], exp_t[i]);
std::fill(exp_t + t, exp_t + t2, 0);
DFT(f, t2);
DFT(exp_t, t2);
for (int i = 0; i != t2; ++i) f[i] = (i64)f[i] * exp_t[i] % mod;
IDFT(f, t2);
std::fill(f + t, f + t2, 0);
}
}
Example
-
Computing
The ordinary approach is polynomial fast exponentiation, with time complexity .
When , we have:
When , let the lowest-degree term of be . Then:
Time complexity .
Polynomial Trigonometric Functions
Given a polynomial , find , and modulo .
Solution
By Euler's formula, , we obtain another expression for trigonometric functions:
Substituting gives:
Writing a program directly from the expressions above gives and modulo . Then can be obtained from .
Code
Polynomial Trigonometric Functions
Notice that we are performing NTT over . Accordingly, the imaginary unit should be replaced by or :
constexpr int MAXN = 262144;
constexpr int mod = 998244353;
constexpr int imgunit = 86583718; /* sqrt(-1) = sqrt(998233452) */
using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;
void polytri(const poly &h, const int n, poly &sin_t, poly &cos_t) {
/* sin(f) = (exp(i * f) - exp(- i * f)) / 2i */
/* cos(f) = (exp(i * f) + exp(- i * f)) / 2 */
/* tan(f) = sin(f) / cos(f) */
assert(h[0] == 0);
static poly_t tri1_t, tri2_t;
for (int i = 0; i != n; ++i) tri2_t[i] = (i64)h[i] * imgunit % mod;
polyexp(tri2_t, n, tri1_t);
polyinv(tri1_t, n, tri2_t);
if (sin_t != nullptr) {
const int invi = fpow(pls(imgunit, imgunit), mod - 2);
for (int i = 0; i != n; ++i)
sin_t[i] = (i64)(tri1_t[i] - tri2_t[i] + mod) * invi % mod;
}
if (cos_t != nullptr) {
for (int i = 0; i != n; ++i) cos_t[i] = div2(pls(tri1_t[i], tri2_t[i]));
}
}
Polynomial Inverse Trigonometric Functions
Given a polynomial , find , and modulo .
Solution
Following the method for computing the polynomial , differentiate and then integrate the inverse trigonometric functions to obtain:
Substituting gives:
The result follows by evaluating these expressions directly.
Code
Polynomial Inverse Trigonometric Functions
constexpr int MAXN = 262144;
constexpr int mod = 998244353;
using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;
void derivative(const poly &h, const int n, poly &f) {
for (int i = 1; i != n; ++i) f[i - 1] = (i64)h[i] * i % mod;
f[n - 1] = 0;
}
void integrate(const poly &h, const int n, poly &f) {
for (int i = n - 1; i; --i) f[i] = (i64)h[i - 1] * inv[i] % mod;
f[0] = 0; /* C */
}
void polyarcsin(const poly &h, const int n, poly &f) {
/* arcsin(f) = ∫ f' / sqrt(1 - f^2) dx */
static poly_t arcsin_t;
const int t = n << 1;
std::copy(h, h + n, arcsin_t);
std::fill(arcsin_t + n, arcsin_t + t, 0);
DFT(arcsin_t, t);
for (int i = 0; i != t; ++i) arcsin_t[i] = sqr(arcsin_t[i]);
IDFT(arcsin_t, t);
arcsin_t[0] = sub(1, arcsin_t[0]);
for (int i = 1; i != n; ++i)
arcsin_t[i] = arcsin_t[i] ? mod - arcsin_t[i] : 0;
polysqrt(arcsin_t, n, f);
polyinv(f, n, arcsin_t);
derivative(h, n, f);
DFT(f, t);
DFT(arcsin_t, t);
for (int i = 0; i != t; ++i) arcsin_t[i] = (i64)f[i] * arcsin_t[i] % mod;
IDFT(arcsin_t, t);
integrate(arcsin_t, n, f);
}
void polyarccos(const poly &h, const int n, poly &f) {
/* arccos(f) = - ∫ f' / sqrt(1 - f^2) dx */
polyarcsin(h, n, f);
for (int i = 0; i != n; ++i) f[i] = f[i] ? mod - f[i] : 0;
}
void polyarctan(const poly &h, const int n, poly &f) {
/* arctan(f) = ∫ f' / (1 + f^2) dx */
static poly_t arctan_t;
const int t = n << 1;
std::copy(h, h + n, arctan_t);
std::fill(arctan_t + n, arctan_t + t, 0);
DFT(arctan_t, t);
for (int i = 0; i != t; ++i) arctan_t[i] = sqr(arctan_t[i]);
IDFT(arctan_t, t);
inc(arctan_t[0], 1);
std::fill(arctan_t + n, arctan_t + t, 0);
polyinv(arctan_t, n, f);
derivative(h, n, arctan_t);
DFT(f, t);
DFT(arctan_t, t);
for (int i = 0; i != t; ++i) arctan_t[i] = (i64)f[i] * arctan_t[i] % mod;
IDFT(arctan_t, t);
integrate(arctan_t, n, f);
}