Polynomial GCD and Resultant

Overview

Polynomial GCD is the backbone of factorization, square-free decomposition, rational-function reduction, and intersection of algebraic curves. Over a field ($\mathbb{Q}[x]$, $\mathrm{GF}(p)[x]$, etc.) the Euclidean algorithm applies directly, but over $\mathbb{Z}[x]$ special care is required to prevent intermediate coefficient explosion.

Related API: gcd(p, q), intPolyGcd, resultant, discriminant, sylvesterMatrix.

Euclidean Algorithm and Coefficient Growth

The Euclidean algorithm over a field reads identically to the integer case:

$$r_0 = f,\; r_1 = g,\; r_{i+1} = r_{i-1} \bmod r_i, \;\gcd(f, g) = r_k \;(r_{k+1} = 0).$$

This works fine over $\mathbb{Q}[x]$, but running a naive pseudo-remainder sequence over $\mathbb{Z}[x]$ causes intermediate coefficients to grow explosively. Knuth's classical example, $f = x^8 + x^6 - 3x^4 - 3x^3 + 8x^2 + 2x - 5$ and $g = 3x^6 + 5x^4 - 4x^2 - 9x + 21$, produces 5-digit intermediate coefficients even though the GCD is the constant 1. In general, intermediate coefficient bit-lengths grow exponentially in the degree.

This is why $\mathbb{Z}[x]$ GCD uses either the subresultant PRS or modular GCD instead of a plain pseudo-remainder sequence.

Subresultant PRS

The subresultant pseudo-remainder sequence (SR-PRS), formalized by Collins (1967) and Brown (1971), computes GCDs over $\mathbb{Z}[x]$ without exponential coefficient growth.

Idea

Pseudo-remainders $r_{i+1} = \mathrm{prem}(r_{i-1}, r_i)$ carry powers of the leading coefficient and are too large. At each step the SR-PRS divides through by an explicit common factor $\beta_i$:

$$r_{i+1} = \frac{\mathrm{prem}(r_{i-1}, r_i)}{\beta_i}.$$

$\beta_i$ has a closed form involving the previous leading coefficient $\ell_i$ and the degree drop (typically $\beta_i = -\ell_{i-1} \cdot \psi_i^{\delta_i}$). Dividing by this factor keeps the result integral while bounding intermediate coefficient bit-lengths polynomially in the input size.

sangi implementation

intPolyGcd(a, b) is implemented with the SR-PRS. pseudoRemainder(f, g) is exposed as the building block (the raw pseudo-remainder). At the final step the content GCD and primitive part are combined to produce the $\mathbb{Z}[x]$ GCD.

Complexity: classical SR-PRS is around $O(n^2 M(n L))$ where $n$ is the degree, $L$ the coefficient bit-length, and $M$ the cost of multi-precision multiplication.

Sylvester Matrix and Resultant

For $f = \sum_{i=0}^n a_i x^i$ and $g = \sum_{i=0}^m b_i x^i$, the Sylvester matrix is the $(m + n) \times (m + n)$ matrix whose upper block stacks $f$'s coefficients shifted across $m$ rows and whose lower block stacks $g$'s coefficients across $n$ rows:

$$\mathrm{Syl}(f, g) = \begin{pmatrix} a_n & a_{n-1} & \cdots & a_0 & & \\ & a_n & a_{n-1} & \cdots & a_0 & \\ & & \ddots & & & \ddots \\ b_m & b_{m-1} & \cdots & b_0 & & \\ & b_m & b_{m-1} & \cdots & b_0 & \\ & & \ddots & & & \ddots \\ \end{pmatrix}.$$

The resultant $\mathrm{Res}(f, g)$ is its determinant:

$$\mathrm{Res}(f, g) = \det \mathrm{Syl}(f, g).$$

Properties

• $\mathrm{Res}(f, g) = 0 \;\Leftrightarrow\; f$ and $g$ share a root in the algebraic closure $\;\Leftrightarrow\; \deg \gcd(f, g) \geq 1$.
• $\mathrm{Res}(f, g) = a_n^m \prod_i g(\alpha_i) = (-1)^{nm} b_m^n \prod_j f(\beta_j)$ where $\alpha_i$ are the roots of $f$ and $\beta_j$ of $g$.
• If $f, g \in \mathbb{Z}[x]$ then $\mathrm{Res}(f, g) \in \mathbb{Z}$ — it is a polynomial in the coefficients only.

Computation

Expanding $\det \mathrm{Syl}$ directly costs $O((n+m)^3)$ and bloats further when coefficients are large. Asymptotically, reading the resultant off the SR-PRS as a by-product is cheaper ($O((n+m)^2 \cdot$ multi-precision coefficient cost$)$): divide the leading coefficient of the last non-zero remainder by the running product of $\beta_i$. sangi's current resultant(f, g) constructs the Sylvester matrix and evaluates its determinant via Gauss elimination ($O((n+m)^3)$); switching to the SR-PRS-fused path is on the roadmap.

Discriminant

The discriminant of $f$ is defined via the resultant with its derivative:

$$\mathrm{Disc}(f) = \frac{(-1)^{n(n-1)/2}}{a_n} \mathrm{Res}(f, f').$$

It generalizes the familiar quadratic discriminant $b^2 - 4ac$.

Properties and uses

• $\mathrm{Disc}(f) = 0 \;\Leftrightarrow\; f$ has a repeated root in its algebraic closure.
• The sign of $\mathrm{Disc}(f)$ can reveal partial parity information about the number of real roots.
• A single resultant computation decides whether $f$ is square-free as a factorization pre-check.

sangi's discriminant(f) internally calls resultant(f, derivative(f)), then applies the sign and division by $a_n$.

// disc(x^2 - 2) = 8
Polynomial<int> f = {-2, 0, 1};
std::cout << discriminant(f) << std::endl;  // 8

Half-GCD (HGCD)

The classical Euclidean / SR-PRS performs $O(n)$ iterations (one per degree drop), each costing $O(n M(n L))$, for a total of $O(n^2 M(n L))$. Half-GCD compresses these into $O(\log n)$ rounds via divide-and-conquer.

Idea

Look only at the high half of the coefficients of $(r_0, r_1)$ and recursively compute a unimodular transformation matrix $\begin{pmatrix} a & b \\ c & d \end{pmatrix}$ that advances the remainder sequence by roughly half the degree. The transformation captures the relation $(r_0, r_1) \to (r_k, r_{k+1})$ with $k \approx n/2$, yielding the GCD and Bezout coefficients in $O(M(n) \log n)$ time.

Complexity: $O(M(n) \log n)$ where $M(n)$ is the cost of polynomial multiplication.

Half-GCD outperforms classical GCD significantly once the degree exceeds a few thousand and FFT-based multiplication becomes effective. sangi's current GCD is SR-PRS-centered; Half-GCD is a planned optional path for larger inputs.

Modular GCD (Brown-Collins)

Compute the GCD over $\mathbb{Z}[x]$ by computing GCDs over $\mathrm{GF}(p)[x]$ for several primes $p_1, p_2, \ldots$ and lifting via the Chinese Remainder Theorem. Due to Brown (1971) and Collins.

Steps

1. Pick a prime $p$ not dividing $\mathrm{lc}(f)$ or $\mathrm{lc}(g)$.
2. Compute $\gcd_p(f, g)$ over $\mathrm{GF}(p)[x]$ via the Euclidean algorithm ($O(n^2)$ in $\mathrm{GF}(p)$).
3. Repeat for several primes, combining via CRT.
4. Stop when the recovered integer coefficients stabilize below the Mignotte bound.

Pros and pitfalls

• No coefficient explosion: every $\mathrm{GF}(p)$ computation stays bounded by $p$.
• Easy to parallelize: the per-prime computations are independent.
• Unlucky primes: a small fraction of primes drop the degree spuriously and must be discarded; the probability is low.

Modular GCD beats SR-PRS when both the coefficient bit-length $L$ and the degree $n$ are large, and is the standard choice in mainstream CAS. sangi's intPolyGcd currently centers on SR-PRS; modular GCD is a planned extension.

Approximate GCD (floating-point)

With double coefficients, exact common factors rarely exist (any small perturbation makes inputs coprime). Approximate GCD extracts an "essential" common factor within a tolerance $\epsilon$.

sangi's approximateGCD(p, q, tol) follows a two-stage strategy:

• Estimate $\deg \gcd$ from the numerical rank of the Sylvester matrix, truncating singular values below $\epsilon$.
• Reconstruct the GCD of that degree by linear least squares.

Signal-processing applications also use methods such as Karmarkar-Lakshman or ApaTools; sangi favors implementation simplicity (rank estimation + LS reconstruction).

Design Decisions

Why SR-PRS as the default

The typical sangi workload (factorization pre-processing, Yun's square-free decomposition, resultant computation) has $n \lesssim 100$ and $L \lesssim 64$, where SR-PRS's small constant outweighs the setup overhead of Half-GCD (FFT) and modular GCD (multiple prime selection, CRT). SR-PRS delivers the best implementation-complexity-versus-performance balance for this range.

Current path and planned SR-PRS switch

Direct expansion of the Sylvester determinant (Gauss elimination) costs $O((n+m)^3)$ and tends to inflate coefficients in $\mathbb{Z}$. Reading the resultant off the SR-PRS as a by-product would fold GCD and resultant computations into a single pass, with a better asymptotic in the multi-precision regime. sangi's current resultant favors implementation simplicity and uses the Sylvester-plus-determinant path; a future release will switch to the SR-PRS-fused route for improved multi-precision performance.

References

• Brown, W. S. (1971). "On Euclid's algorithm and the computation of polynomial greatest common divisors". Journal of the ACM, 18(4), 478–504.
• Collins, G. E. (1967). "Subresultants and reduced polynomial remainder sequences". Journal of the ACM, 14(1), 128–142.
• Knuth, D. E. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley, 1997.
• von zur Gathen, J. and Gerhard, J. Modern Computer Algebra. 3rd ed. Cambridge University Press, 2013.