viii. Polynomial GCD and Resultant

The Euclidean algorithm for polynomials over a field: gcd, Bézout coefficients and inverses modulo a polynomial in $O(n{\log }^{2}n)$ by half-gcd, and the resultant $\mathrm{Res}(A,B)$ (zero iff $A$, $B$ share a root) with its uses: eliminating a variable, discriminants, circulant determinants. Uses the core of NTT and Power Series (polyMul, polyDivMod).

Task	Method	Cost
$\gcd (A,B)$, Bézout $UA+VB=\gcd$	schoolbook Euclid	$O(nm)$
same, large degrees	half-gcd (Knuth–Schönhage)	$O(n{\log }^{2}n)$
$A^{-1}\mathrm{mod}B$	Bézout coefficient $U$ when $\gcd =1$ (Library Checker “Inv of Polynomials”)	$O(n{\log }^{2}n)$
$\mathrm{Res}(A,B)$	Euclid with leading-coefficient bookkeeping	$O(nm)$ (half-gcd version $O(n{\log }^{2}n)$)
discriminant $\mathrm{disc}(f)=\frac{(-1{)}^{n(n-1)/2}}{\mathrm{lc}(f)}\mathrm{Res}(f,f^{\prime })$	resultant	$O(n^2)$
squarefree part $f/\gcd (f,f^{\prime })$	gcd	as above
eliminate $y$ from $P(x,y)=Q(x,y)=0$	${\mathrm{Res}}_y(P,Q)$ is a polynomial in $x$ vanishing at the common solutions	interpolate over $x$

$\mathrm{Res}(A,B)=\mathrm{lc}(A{)}^{\mathrm{deg}B}{\prod }_{A(\alpha )=0}B(\alpha )=(-1{)}^{\mathrm{deg}A\mathrm{deg}B}\mathrm{lc}(B{)}^{\mathrm{deg}A}{\prod }_{B(\beta )=0}A(\beta )=\det (\text{Sylvester matrix}).$

Half-gcd. One Euclid step maps $(a,b)\mapsto (b,a-qb)$, a $2\times 2$ matrix $\left(\begin{array}{cc}0& 1\\ 1& -q\end{array}\right)$. The quotients of the first steps depend only on the top coefficients: from the top halves of $a$ and $b$ one computes recursively the product $M$ of all steps that bring the degree below $\mathrm{deg}a/2$, applies it, does one step, and recurses on the top half again; the full gcd repeats half-gcd.

struct PolyEuclid {                                         // 2x2 polynomial matrices for Euclid steps
    using Pair = array<Poly, 2>; using M = array<Pair, 2>;
    static Poly strip(Poly a) { while (!a.empty() && a.back() == mint(0)) a.pop_back(); return a; }
    static Poly mul(const Poly& a, const Poly& b) { return a.empty() || b.empty() ? Poly() : polyMul(a, b); }
    static Poly add(Poly a, const Poly& b, bool neg = false) { if (sz(a) < sz(b)) a.resize(sz(b)); for (int i = 0; i < sz(b); i++) a[i] += neg ? -b[i] : b[i]; return strip(a); }
    static M identity() { return {Pair{Poly{1}, Poly{}}, Pair{Poly{}, Poly{1}}}; }
    static Pair apply(const M& m, const Pair& v) { Pair r; for (int i = 0; i < 2; i++) r[i] = add(mul(m[i][0], v[0]), mul(m[i][1], v[1])); return r; }
    static M compose(const M& a, const M& b) { M r; for (int i = 0; i < 2; i++) for (int k = 0; k < 2; k++) r[i][k] = add(mul(a[i][0], b[0][k]), mul(a[i][1], b[1][k])); return r; }
    static Pair cut(Pair a, int k) { for (auto& p : a) p = sz(p) <= k ? Poly() : Poly(p.begin() + k, p.end()); return a; }   // divide by x^k
    static void step(M& m, Pair& a) {                        // (a0, a1) -> (a1, a0 - q a1)
        Poly q = strip(polyDivMod(a[0], a[1]).first);
        Poly r = add(a[0], mul(q, a[1]), true);
        a = {a[1], r};
        m = {m[1], Pair{add(m[0][0], mul(q, m[1][0]), true), add(m[0][1], mul(q, m[1][1]), true)}};
    }
    static M halfGcd(Pair a) {                              // deg a0 > deg a1 >= 0
        int m = sz(a[0]) / 2;
        if (sz(a[1]) <= m) return identity();
        M R = halfGcd(cut(a, m)); a = apply(R, a);
        if (sz(a[1]) <= m) return R;
        step(R, a);
        if (sz(a[1]) <= m) return R;
        int l = sz(a[0]) - 1, k = 2 * m - l;
        return compose(halfGcd(cut(a, k)), R);
    }
    static M euclid(Pair a) {                               // all steps until a1 = 0
        M m0 = halfGcd(a); a = apply(m0, a);
        if (a[1].empty()) return m0;
        step(m0, a);
        if (a[1].empty()) return m0;
        return compose(euclid(a), m0);
    }
    static M run(const Poly& a, const Poly& b) {            // M (a, b) = (gcd, 0) up to a constant
        M res = identity(); Pair v = {strip(a), strip(b)};
        if (v[1].empty()) return res;
        step(res, v);
        if (!v[1].empty()) res = compose(euclid(v), res);
        return res;
    }
};
Poly polyGcd(Poly a, Poly b) {                              // monic gcd, {} if both zero
    a = PolyEuclid::strip(a), b = PolyEuclid::strip(b);
    if (b.empty()) swap(a, b);
    if (b.empty()) return {};
    auto m = PolyEuclid::run(a, b);
    Poly g = PolyEuclid::add(PolyEuclid::mul(m[0][0], a), PolyEuclid::mul(m[0][1], b));
    mint inv = g.back().inv(); for (auto& x : g) x *= inv;
    return g;
}
Poly polyInverseMod(Poly a, Poly b) {                       // x with a x = 1 mod b, deg x < deg b; {} if gcd != 1
    a = PolyEuclid::strip(a), b = PolyEuclid::strip(b);
    if (b.empty() || a.empty()) return {};
    auto m = PolyEuclid::run(a, b);
    Poly g = PolyEuclid::add(PolyEuclid::mul(m[0][0], a), PolyEuclid::mul(m[0][1], b));
    if (sz(g) != 1) return {};
    Poly x = m[0][0]; mint inv = g[0].inv(); for (auto& c : x) c *= inv;
    if (sz(x) >= sz(b)) x = polyDivMod(x, b).second;
    return x;
}
 
mint resultant(Poly a, Poly b) {                           // O(deg a deg b), any lengths
    a = PolyEuclid::strip(a), b = PolyEuclid::strip(b);
    if (a.empty() || b.empty()) return 0;
    mint res = 1;
    while (true) {
        if (sz(b) == 1) { for (int i = 0; i + 1 < sz(a); i++) res *= b[0]; return res; }   // Res(A, c) = c^deg A
        if (sz(a) == 1) { for (int i = 0; i + 1 < sz(b); i++) res *= a[0]; return res; }   // Res(c, B) = c^deg B
        Poly c = polyDivMod(a, b).second;
        if (c.empty()) return 0;
        int da = sz(a) - 1, db = sz(b) - 1, dc = sz(c) - 1;
        if ((ll)da * db % 2) res = -res;
        for (int i = 0; i < da - dc; i++) res *= b.back();
        a = b, b = c;
    }
}
mint discriminant(const Poly& f) {                          // (-1)^{n(n-1)/2} Res(f, f') / lc(f), n = deg f >= 1
    Poly g = PolyEuclid::strip(f); int n = sz(g) - 1;
    mint r = resultant(g, derivative(g)) / g.back();
    return ((ll)n * (n - 1) / 2) % 2 ? -r : r;
}

Proof: half-gcd, Bézout and the resultant recursion

Euclid as matrices. Each step multiplies the pair by $\left(\begin{array}{cc}0& 1\\ 1& -q\end{array}\right)$ (invertible, determinant $-1$), so after all steps $M(a,b{)}^T=(g,0{)}^T$ with $g$ a gcd and the first row of $M$ the Bézout coefficients; step updates $M$ by left multiplication. Half-gcd. Write $a_i=h_i{x}^m+{\ell }_i$ with $\mathrm{deg}{\ell }_i<m$. Key lemma (von zur Gathen–Gerhard, Modern Computer Algebra, §11.1): the Euclidean quotients of $(h_0,h_1)$, as long as the remainders keep degree at least $\frac{1}{2}\mathrm{deg}{h}_0$, coincide with the corresponding quotients of $(a_0,a_1)$, because a quotient $q$ with $\mathrm{deg}q=d$ is determined by the top $d+1$ coefficients of dividend and divisor, and the low parts ${\ell }_i$ only influence coefficients below $x^m$ times the growth of the previous quotients. Hence the matrix $R$ returned for $\mathrm{cut}(a,m)$ consists of genuine Euclid steps of $a$; after applying it and one explicit step, the second recursive call on the top part (cut at $k=2m-l$) again performs only genuine steps, and the procedure stops exactly when the remainder degree drops to at most $m$ (the sz <= m tests). Each call makes two recursive calls on inputs of half the size plus $O(n\log n)$ multiplication work, so $T(n)=2T(n/2)+O(n\log n)=O(n{\log }^{2}n)$. euclid alternates half-gcd and single steps until the remainder is 0 (the degree halves each round, so the total is still $O(n{\log }^{2}n)$). Inverse: if $Ua+Vb=c\ne 0$ constant then $(U/c)a\equiv 1(\mathrm{mod}b)$, reduced mod $b$.

Resultant. From the product formula $\mathrm{Res}(A,B)=\mathrm{lc}(A{)}^{\mathrm{deg}B}{\prod }_{A(\alpha )=0}B(\alpha )$ (roots in a splitting field): $\mathrm{Res}(A,B)=(-1{)}^{\mathrm{deg}A\mathrm{deg}B}\mathrm{Res}(B,A)$, and for $A=QB+C$, $\mathrm{Res}(B,A)=\mathrm{lc}(B{)}^{\mathrm{deg}A}{\prod }_{B(\beta )=0}C(\beta )=\mathrm{lc}(B{)}^{\mathrm{deg}A-\mathrm{deg}C}\mathrm{Res}(B,C)$; a constant $c$ gives $\mathrm{Res}(A,c)=c^{\mathrm{deg}A}$; a zero remainder with $\mathrm{deg}B\ge 1$ means a common factor, so the resultant is 0. The product formula equals the Sylvester determinant (both are the unique polynomial in the coefficients that vanishes exactly on common roots, with the same normalisation; or directly by the standard row-reduction argument). Discriminant: $\mathrm{Res}(f,f^{\prime })={\mathrm{lc}}^{n-1}{\prod }_i{f}^{\prime }({\alpha }_i)$ and $f^{\prime }({\alpha }_i)=\mathrm{lc}{\prod }_{j\ne i}({\alpha }_i-{\alpha }_j)$, so $\mathrm{Res}(f,f^{\prime })={\mathrm{lc}}^{2n-1}{\prod }_{i\ne j}({\alpha }_i-{\alpha }_j)=(-1{)}^{n(n-1)/2}\mathrm{lc}\cdot \mathrm{disc}(f)$ with $\mathrm{disc}={\mathrm{lc}}^{2n-2}{\prod }_{i<j}({\alpha }_i-{\alpha }_j{)}^2$. Checked: half-gcd against schoolbook Euclid for random and structured (common factor) inputs; resultant against the Sylvester determinant; discriminant of quadratics $b^2-4ac$.

Related notes

• Structured Matrices (circulant determinant via resultant)
• Polynomial Factorization (gcd with $x^{p^d}-x$)
• GCD and Extended Euclid (the integer analogue)