xi. Polynomial Factorization

Factor a polynomial into irreducibles. Over a finite field this is polynomial time (randomised Cantor–Zassenhaus, or deterministic Berlekamp for small fields); over $\mathbb{Z}$ it reduces to the finite-field case plus Hensel lifting and recombination (LLL makes it polynomial). The contest use is roots modulo $p$. Uses Poly, polyMul, polyDivMod, derivative (NTT, Power Series) and polyGcd (GCD and Resultant), over ${\mathbb{F}}_p$ with $p=998244353$.

The pipeline over ${\mathbb{F}}_q$: (1) square-free factorisation $f={\prod }_i{g}_i^i$ with $g_i$ square-free, from gcds with $f^{\prime }$; (2) distinct-degree factorisation: $\gcd (f,x^{q^d}-x)$ collects exactly the irreducible factors of degree $d$ once smaller degrees are removed; (3) equal-degree factorisation: split a product of irreducibles of the same degree $d$ randomly with $\gcd (g,h^{(q^d-1)/2}-1)$ (for $q=2^m$ use the trace $h+h^2+\cdots +h^{2^{dm-1}}$).

Algorithm	Type	Time	Best for
Berlekamp	deterministic	$O(n^3+q{n}^2)$	small $q$ (e.g. ${\mathbb{F}}_2$); gives the number of factors $n-\mathrm{rank}(Q-I)$ first
Cantor–Zassenhaus	Las Vegas	$\tilde{O}(n^2\log q)$	large $q$: the standard
von zur Gathen–Shoup	randomised	$\tilde{O}(n^{1.5}\log q)$ (baby-step giant-step DDF)	asymptotically faster DDF
Kaltofen–Shoup	randomised	$\tilde{O}(n^{1.815}\log q)$	theoretical
over $\mathbb{Z}$	mod $p$ + Hensel lift past twice the Mignotte bound + recombination (LLL: Lenstra–Lenstra–Lovász)	polynomial with LLL; exponential recombination in bad cases (Swinnerton-Dyer)

Test / count	Method
irreducible over ${\mathbb{F}}_p$ (Rabin)	$x^{p^n}\equiv x(\mathrm{mod}f)$ and $\gcd (f,x^{p^{n/r}}-x)=1$ for every prime $r\mid n$
irreducible via Berlekamp	$\mathrm{rank}(Q-I)=n-1$
over $\mathbb{Z}$: Eisenstein	a prime $p$ divides all coefficients but the leading one and $p^2\nmid a_0$
over $\mathbb{Z}$	irreducible mod some $p\nmid \mathrm{lc}$ (with the same degree) $\Rightarrow$ irreducible over $\mathbb{Q}$ (converse fails)
cyclotomic ${\Phi }_n$	always irreducible over $\mathbb{Q}$
squarefree	$\gcd (f,f^{\prime })=1$
monic irreducibles of degree $n$ over ${\mathbb{F}}_q$	$\frac{1}{n}{\sum }_{d\mid n}\mu (d)q^{n/d}$ (the Lyndon-word formula)

mt19937_64 frng(7);
Poly polyModP(const Poly& a, const Poly& f) { return polyDivMod(a, f).second; }
Poly mulModP(const Poly& a, const Poly& b, const Poly& f) { return polyModP(polyMul(a, b), f); }
Poly powModP(Poly b, ll e, const Poly& f) { Poly r = {1}; b = polyModP(b, f); for (; e; e >>= 1, b = mulModP(b, b, f)) if (e & 1) r = mulModP(r, b, f); return polyModP(r, f); }
Poly polySub1(Poly t) { if (t.empty()) t = {0}; t[0] -= 1; while (!t.empty() && t.back() == mint(0)) t.pop_back(); return t; }
const ll Q = 998244353;
 
void equalDegreeSplit(const Poly& g, int d, vector<Poly>& out) {   // g monic squarefree, all factors of degree d
    int n = sz(g) - 1; if (n == 0) return;
    if (n == d) { out.push_back(g); return; }
    while (true) {
        Poly h(n); for (auto& x : h) x = frng() % Q;
        Poly cur = polyModP(h, g), prod = cur;              // h^{(q^d - 1)/2} = (h h^q ... h^{q^{d-1}})^{(q-1)/2}
        for (int i = 1; i < d; i++) cur = powModP(cur, Q, g), prod = mulModP(prod, cur, g);
        Poly r = polyGcd(g, polySub1(powModP(prod, (Q - 1) / 2, g)));
        int dr = sz(r) - 1;
        if (dr > 0 && dr < n) { equalDegreeSplit(r, d, out); equalDegreeSplit(polyDivMod(g, r).first, d, out); return; }
    }
}
vector<pair<Poly, int>> squareFreeDecomposition(Poly f) {  // Yun, monic f, deg f < p: f = prod g_i^i
    vector<pair<Poly, int>> res; Poly df = derivative(f), a = polyGcd(f, df);
    if (sz(a) == 0) a = {1};
    Poly b = polyDivMod(f, a).first, c = polyDivMod(df, a).first;
    Poly d = c; { Poly db = derivative(b); if (sz(d) < sz(db)) d.resize(sz(db)); for (int i = 0; i < sz(db); i++) d[i] -= db[i]; }
    for (int i = 1; sz(b) > 1; i++) {
        Poly ai = polyGcd(b, d); if (ai.empty()) ai = b;
        if (sz(ai) > 1) res.push_back({ai, i});
        b = polyDivMod(b, ai).first; c = polyDivMod(d, ai).first;
        Poly db = derivative(b); d = c; if (sz(d) < sz(db)) d.resize(sz(db)); for (int j = 0; j < sz(db); j++) d[j] -= db[j];
        while (!d.empty() && d.back() == mint(0)) d.pop_back();
    }
    return res;
}
vector<pair<Poly, int>> factorModP(Poly f) {              // monic irreducible factors with multiplicities (deg f < p)
    while (!f.empty() && f.back() == mint(0)) f.pop_back();
    vector<pair<Poly, int>> res; if (sz(f) <= 1) return res;
    mint lc = f.back().inv(); for (auto& x : f) x *= lc;
    for (auto [g, mult] : squareFreeDecomposition(f)) {
        Poly rest = g, xq = {0, 1};
        for (int d = 1; 2 * d <= sz(rest) - 1; d++) {        // distinct-degree factorisation
            xq = powModP(xq, Q, rest);                      // x^{q^d} mod rest
            Poly part = polyGcd(rest, [&] { Poly t = xq; if (sz(t) < 2) t.resize(2); t[1] -= 1; while (!t.empty() && t.back() == mint(0)) t.pop_back(); return t; }());
            if (sz(part) > 1) {
                vector<Poly> pieces; equalDegreeSplit(part, d, pieces);
                for (auto& p : pieces) res.push_back({p, mult});
                rest = polyDivMod(rest, part).first; xq = polyModP(xq, rest);
            }
        }
        if (sz(rest) > 1) res.push_back({rest, mult});       // remaining factor is irreducible
    }
    return res;
}
vector<mint> rootsModP(Poly f) {                            // distinct roots of f in F_p
    while (!f.empty() && f.back() == mint(0)) f.pop_back();
    if (sz(f) <= 1) return {};
    Poly xp = powModP({0, 1}, Q, f); if (sz(xp) < 2) xp.resize(2); xp[1] -= 1;   // x^p - x mod f
    while (!xp.empty() && xp.back() == mint(0)) xp.pop_back();
    Poly g = xp.empty() ? f : polyGcd(f, xp);                // product of the distinct linear factors
    { mint lc = g.back().inv(); for (auto& x : g) x *= lc; }
    vector<Poly> lin; equalDegreeSplit(g, 1, lin);
    vector<mint> res; for (auto& l : lin) res.push_back(-l[0]);
    sort(all(res), [](mint a, mint b) { return a.v < b.v; }); return res;
}
bool isIrreducible(Poly f) {                                // Rabin's test
    while (!f.empty() && f.back() == mint(0)) f.pop_back();
    int n = sz(f) - 1; if (n <= 0) return false;
    mint lc = f.back().inv(); for (auto& x : f) x *= lc;
    vector<int> primes; for (int m = n, r = 2; m > 1; r++) if (m % r == 0) { primes.push_back(r); while (m % r == 0) m /= r; }
    auto xPowQPow = [&](int k) { Poly t = {0, 1}; for (int i = 0; i < k; i++) t = powModP(t, Q, f); return t; };   // x^{q^k} mod f
    Poly full = xPowQPow(n); Poly xx = polyModP({0, 1}, f);
    if (full != xx) return false;
    for (int r : primes) { Poly t = xPowQPow(n / r); if (sz(t) < 2) t.resize(2); t[1] -= 1; while (!t.empty() && t.back() == mint(0)) t.pop_back(); if (t.empty() || sz(polyGcd(f, t)) > 1) return false; }
    return true;
}

Proof: the finite-field facts behind SFF, DDF, EDF, root finding and Rabin

$x^{q^d}-x$. The roots of $x^{q^d}-x$ in an algebraic closure are exactly the $q^d$ elements of ${\mathbb{F}}_{q^d}$, each simple (derivative $-1$). A monic irreducible $\pi$ of degree $e$ has its roots in ${\mathbb{F}}_{q^e}$, which lies inside ${\mathbb{F}}_{q^d}$ iff $e\mid d$. So $x^{q^d}-x={\prod }_{\mathrm{deg}\pi \mid d}\pi$. DDF: after the factors of degree $<d$ are divided out of a square-free $f$, $\gcd (f,x^{q^d}-x)$ is the product of its degree-$d$ factors; once $2d>\mathrm{deg}(\text{rest})$, the rest has no factor of degree $\le d$ and is irreducible. SFF (Yun): for $f=\prod g_i^i$ with pairwise coprime square-free $g_i$ in characteristic $p>\mathrm{deg}f$, $f^{\prime }={\sum }_{i}i{g}_i^{i-1}{g}_i^{\prime }{\prod }_{j\ne i}{g}_j^j$, so $\gcd (f,f^{\prime })=\prod g_i^{i-1}$; Yun’s iteration maintains $b={\prod }_{j\ge i}{g}_j$ and $d=c-b^{\prime }=g_i\cdot (\dots )$ with $\gcd (b,d)=g_i$, peeling one multiplicity per round. (In characteristic $p\le \mathrm{deg}f$, $f^{\prime }=0$ means $f$ is a $p$-th power: take the $p$-th root of the coefficients.)

EDF (Cantor–Zassenhaus). By CRT, ${\mathbb{F}}_q[x]/(g)\cong {\prod }_{i=1}^k{\mathbb{F}}_{q^d}$. For uniform random $h$ the components are independent uniform elements of ${\mathbb{F}}_{q^d}$, and $u\mapsto u^{(q^d-1)/2}$ maps ${\mathbb{F}}_{q^d}^{\ast }$ onto $\{\pm 1\}$, each with probability $\frac{1}{2}$ (half of the cyclic group are squares). The gcd with $h^{(q^d-1)/2}-1$ is the product of the factors whose component is $+1$, a non-trivial split unless all components agree, which has probability at most $2^{1-k}+O(k/q^d)\le \frac{1}{2}+o(1)$ for $k\ge 2$. The exponent is computed as $(h{h}^q\cdots h^{q^{d-1}}{)}^{(q-1)/2}$ because $\frac{q^d-1}{2}=\frac{q-1}{2}(1+q+\cdots +q^{d-1})$. Roots: $\gcd (f,x^p-x)$ is the product of the distinct linear factors of $f$ ($d=1$ above), which EDF splits; this generalises Tonelli–Shanks ($f=x^2-a$). Rabin: $f$ of degree $n$ is irreducible iff it divides $x^{q^n}-x$ (all factors have degree dividing $n$) and has no factor of degree dividing $n/r$ for a prime $r\mid n$ (every proper divisor of $n$ divides some $n/r$). All functions were checked: factor products reproduce $f$, every factor passes Rabin and has the degree reported, roots match polynomials built from known roots, and Rabin agrees with brute-force trial division over ${\mathbb{F}}_5$-sized analogues.

Berlekamp

For square-free $f=f_1\cdots f_k$, CRT gives ${\mathbb{F}}_q[x]/(f)\cong {\prod }_i{\mathbb{F}}_q[x]/(f_i)$. The map $\varphi (g)=g^q-g\mathrm{mod}f$ is ${\mathbb{F}}_q$-linear and its kernel, the Berlekamp subalgebra, consists exactly of the elements that are constant in every CRT component, so it has dimension $k$. Build the matrix $Q$ whose row $i$ is $x^{qi}\mathrm{mod}f$; the kernel of $Q-I$ gives $k=n-\mathrm{rank}(Q-I)$ and a basis $1=g_1,g_2,\dots ,g_k$. For a non-constant $g$ in the kernel and $c\in {\mathbb{F}}_q$, $\gcd (f,g-c)$ is the product of the $f_i$ whose component of $g$ is $c$: a non-trivial factor for some $c$. Building $Q$ is $O(n^2\log q)$, the kernel $O(n^3)$, splitting $O(q{n}^2)$: linear in $q$, fine for small primes (coding theory over ${\mathbb{F}}_2$), hopeless for large ones, where Cantor–Zassenhaus-style random splitting replaces it. Berlekamp designed it for error-correcting codes, like his Berlekamp–Massey algorithm (Linear Recurrences).

using PolyP = vector<ll>;                                  // coefficients mod a small prime p, low degree first
PolyP trimP(PolyP a) { while (!a.empty() && a.back() == 0) a.pop_back(); return a; }
PolyP modP(PolyP a, const PolyP& b, ll p) {                // b monic
    a = trimP(a); int m = sz(b) - 1;
    for (int i = sz(a) - 1; i >= m; i--) { ll c = a[i]; if (!c) continue; for (int j = 0; j <= m; j++) a[i - m + j] = ((a[i - m + j] - c * b[j]) % p + p) % p; }
    if (sz(a) > m) a.resize(m); return trimP(a);
}
ll invP(ll a, ll p) { ll r = 1, e = p - 2; a %= p; for (; e; e >>= 1, a = a * a % p) if (e & 1) r = r * a % p; return r; }
PolyP monicP(PolyP a, ll p) { a = trimP(a); if (a.empty()) return a; ll iv = invP(a.back(), p); for (auto& x : a) x = x * iv % p; return a; }
PolyP gcdP(PolyP a, PolyP b, ll p) { a = trimP(a), b = trimP(b); while (!b.empty()) { b = monicP(b, p); a = modP(a, b, p); swap(a, b); } return monicP(a, p); }
PolyP mulP(const PolyP& a, const PolyP& b, ll p) { if (a.empty() || b.empty()) return {}; PolyP c(sz(a) + sz(b) - 1, 0); for (int i = 0; i < sz(a); i++) for (int j = 0; j < sz(b); j++) c[i + j] = (c[i + j] + a[i] * b[j]) % p; return c; }
PolyP divP(PolyP a, const PolyP& b, ll p) {                // exact quotient, b monic
    a = trimP(a); int m = sz(b) - 1; if (sz(a) - 1 < m) return {};
    PolyP q(sz(a) - m, 0);
    for (int i = sz(a) - 1; i >= m; i--) { ll c = a[i]; q[i - m] = c; if (!c) continue; for (int j = 0; j <= m; j++) a[i - m + j] = ((a[i - m + j] - c * b[j]) % p + p) % p; }
    return q;
}
vector<PolyP> berlekampFactor(PolyP f, ll p) {             // f monic square-free over F_p (small p); irreducible factors
    int n = sz(f) - 1; if (n <= 1) return {f};
    vector<vector<ll>> M(n, vector<ll>(n, 0));             // row i: x^{p i} mod f, minus identity
    PolyP xp = {1}, base = modP([&] { PolyP t(p + 1, 0); t[p] = 1; return t; }(), f, p);
    for (int i = 0; i < n; i++) { for (int j = 0; j < sz(xp); j++) M[i][j] = xp[j]; M[i][i] = (M[i][i] + p - 1) % p; xp = modP(mulP(xp, base, p), f, p); }
    // kernel of M^T: vectors v with sum_i v_i (row_i) = 0, i.e. g = sum v_i x^i with g^p = g mod f
    vector<vector<ll>> A(n, vector<ll>(n)); for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) A[j][i] = M[i][j];
    vector<int> where(n, -1); int row = 0;
    for (int col = 0; col < n && row < n; col++) {
        int piv = -1; for (int r = row; r < n; r++) if (A[r][col]) { piv = r; break; }
        if (piv < 0) continue;
        swap(A[piv], A[row]); ll iv = invP(A[row][col], p); for (auto& x : A[row]) x = x * iv % p;
        for (int r = 0; r < n; r++) if (r != row && A[r][col]) { ll c = A[r][col]; for (int j = 0; j < n; j++) A[r][j] = ((A[r][j] - c * A[row][j]) % p + p) % p; }
        where[col] = row++;
    }
    vector<PolyP> basis;
    for (int fcol = 0; fcol < n; fcol++) if (where[fcol] == -1) {
        PolyP v(n, 0); v[fcol] = 1; for (int j = 0; j < n; j++) if (where[j] != -1) v[j] = (p - A[where[j]][fcol]) % p;
        basis.push_back(trimP(v));
    }
    int k = sz(basis); vector<PolyP> fs = {f};
    for (auto& g : basis) {
        if (sz(fs) == k) break;
        if (sz(g) <= 1) continue;                           // constants split nothing
        vector<PolyP> nxt;
        for (auto& h : fs) {
            if (sz(h) <= 2) { nxt.push_back(h); continue; }
            for (ll c = 0; c < p && sz(h) > 1; c++) {
                PolyP gc = g; gc[0] = (gc[0] - c % p + p) % p; PolyP d = gcdP(h, gc, p);
                if (sz(d) > 1 && sz(d) < sz(h)) { nxt.push_back(d); h = divP(h, d, p); }
            }
            if (sz(h) > 1) nxt.push_back(h);
        }
        fs = nxt;
    }
    return fs;
}
ll countMonicIrreducible(int n, ll q) {                    // (1/n) sum_{d | n} mu(d) q^{n/d}, exact for small values
    ll s = 0;
    for (int d = 1; d <= n; d++) if (n % d == 0) {
        int mu = 1, x = d; for (int r = 2; r <= x; r++) if (x % r == 0) { x /= r; if (x % r == 0) { mu = 0; break; } mu = -mu; }
        ll pw = 1; for (int i = 0; i < n / d; i++) pw *= q; s += mu * pw;
    }
    return s / n;
}

Proof: the Berlekamp subalgebra, the splitting step and the irreducible count

Kernel. In ${\prod }_i{\mathbb{F}}_q[x]/(f_i)\cong {\prod }_i{\mathbb{F}}_{q^{d_i}}$, $g^q=g$ holds iff every component satisfies $u^q=u$, i.e. lies in the prime-power subfield ${\mathbb{F}}_q$ (the roots of $u^q-u$). So the kernel is ${\mathbb{F}}_q^k$, of dimension $k$. Since $g\mapsto g^q$ is linear over ${\mathbb{F}}_q$ (Frobenius), it is the matrix $Q$ in the basis $1,x,\dots ,x^{n-1}$, and the kernel of $Q-I$ (acting on coefficient row vectors) is computed by elimination on the transpose. Splitting. For $g$ in the kernel with components $(c_1,\dots ,c_k)$, $f_i\mid g-c$ iff $c_i=c$, so $\gcd (f,g-c)={\prod }_{c_i=c}{f}_i$. Two factors $f_i\ne f_j$ are separated by some basis vector (otherwise every kernel element would have $c_i=c_j$, contradicting $\dim =k$), so processing all basis vectors with all $c$ separates every pair, and the loop stops when $k$ factors are found. Counting irreducibles: comparing degrees in $x^{q^n}-x={\prod }_{\mathrm{deg}\pi \mid n}\pi$ gives $q^n={\sum }_{d\mid n}d{I}_d$, and Möbius inversion gives $I_n=\frac{1}{n}{\sum }_{d\mid n}\mu (d)q^{n/d}$. Checked over ${\mathbb{F}}_2,{\mathbb{F}}_3,{\mathbb{F}}_5$ against brute-force factorisation by trial division and against enumerated irreducible counts.

Cyclotomic polynomials

$x^n-1={\prod }_{d\mid n}{\Phi }_d(x)$ is the factorisation over $\mathbb{Q}$; ${\Phi }_n$ (degree $\phi (n)$) has the primitive $n$-th roots of unity as roots, integer coefficients, and is irreducible over $\mathbb{Q}$. ${\Phi }_n(x)={\prod }_{d\mid n}(x^d-1{)}^{\mu (n/d)}$. Coefficients are in $\{-1,0,1\}$ for $n<105$ (${\Phi }_{105}$ has a $-2$) but grow for $n$ with many odd prime factors. ${\Phi }_p=1+x+\cdots +x^{p-1}$; ${\Phi }_n(1)=p$ if $n=p^k$, else 1 ($n>1$).

vector<ll> cyclotomic(int n) {                              // coefficients of Phi_n, O(n d(n))
    vector<int> mu(n + 1, 1); vector<char> comp(n + 1, 0);
    for (int p = 2; p <= n; p++) if (!comp[p]) for (int j = p; j <= n; j += p) { if (j > p) comp[j] = 1; mu[j] = (j / p) % p == 0 ? 0 : -mu[j]; }
    vector<ll> poly = {1};
    auto mulXd = [&](int d) { vector<ll> r(sz(poly) + d, 0); for (int i = 0; i < sz(poly); i++) r[i + d] += poly[i], r[i] -= poly[i]; poly = r; };   // * (x^d - 1)
    auto divXd = [&](int d) { int m = sz(poly); vector<ll> r(m - d, 0); for (int i = m - 1; i >= d; i--) { r[i - d] = poly[i]; poly[i - d] += poly[i]; } poly = r; };   // / (x^d - 1), exact
    for (int d = 1; d <= n; d++) if (n % d == 0 && mu[n / d] == 1) mulXd(d);
    for (int d = 1; d <= n; d++) if (n % d == 0 && mu[n / d] == -1) divXd(d);
    return poly;
}
vector<pair<int, vector<ll>>> factorXnMinus1(int n) { vector<pair<int, vector<ll>>> r; for (int d = 1; d <= n; d++) if (n % d == 0) r.push_back({d, cyclotomic(d)}); return r; }

Proof: the cyclotomic product formula

Every $n$-th root of unity is a primitive $d$-th root for exactly one $d\mid n$, so $x^n-1={\prod }_{d\mid n}{\Phi }_d$. Taking logarithms (or working in the multiplicative group of non-zero rational functions) and applying Möbius inversion gives ${\Phi }_n={\prod }_{d\mid n}(x^d-1{)}^{\mu (n/d)}$. Multiplying the factors with $\mu =1$ first and then dividing by those with $\mu =-1$ keeps everything a polynomial (each partial quotient is a product of cyclotomic polynomials) and the divisions by the monic $x^d-1$ are exact with integer coefficients (the division loop is synthetic division by $x^d-1$). Irreducibility over $\mathbb{Q}$ is Gauss’s / Dedekind’s theorem (quoted). Verified: ${\prod }_{d\mid n}{\Phi }_d=x^n-1$ for $n\le 300$, $\mathrm{deg}{\Phi }_n=\phi (n)$, and the $-2$ in ${\Phi }_{105}$.

Related notes

• Modular Square Roots (Tonelli–Shanks is $f=x^2-a$)
• GCD and Resultant
• Lyndon Words (the same counting formula)