8.4. The half-gcd algorithm🔗

Euclid's algorithm divides r₀ = a by r₁ = b, then r₁ by the remainder r₂, and so on, with rᵢ₊₁ = rᵢ₋₁ - qᵢ rᵢ. Each step is a multiplication by the matrix [[0, 1], [1, -qᵢ]], so the pair (rᵢ, rᵢ₊₁) is the product of the matrices so far applied to (a, b). When a and b have degree about n, the quotients in the first half of the algorithm, until the degree has dropped by n/2, depend only on the top halves of the coefficients of a and b. The half-gcd algorithm computes the product of those quotient matrices recursively from the top halves, applies it to a and b, and repeats the process on the result. With fast multiplication, this costs O(M(n) log n), where M(n) is the cost of one product.

This implementation does not trust the top halves. It computes the matrix from the top halves, applies it to the whole of a and b, and accepts it only if the two resulting polynomials have the degrees that the top halves predict. Otherwise it discards the matrix and takes one ordinary Euclidean step with Hex.DensePoly.divModWith. Accepted matrices are proved to be products of Euclid's own quotient matrices, so the algorithm returns exactly the gcd and Bézout coefficients of HexPoly's xgcd, which are not normalized to make the gcd monic. The gcd functions, interpolation and Padé approximation need a field of coefficients. Matrices are values of type Hex.DensePoly.GcdStep.

{ gcd := #p[-1, 1], left := #p[1], right := #p[0, -1] }

Here (x³ - 1) · 1 + (x² - 1) · (-x) = x - 1.

Hex.DensePoly.xgcdWith_eq.{u} {F : Type u} [DecidableEq F]
  [Lean.Grind.Field F] (plan : Hex.DensePoly.MulPlan F)
  (p q : Hex.DensePoly F) : Hex.DensePoly.xgcdWith plan p q = p.xgcd q
Hex.DensePoly.xgcdWith_eq.{u} {F : Type u}
  [DecidableEq F] [Lean.Grind.Field F]
  (plan : Hex.DensePoly.MulPlan F)
  (p q : Hex.DensePoly F) :
  Hex.DensePoly.xgcdWith plan p q =
    p.xgcd q

Hex.DensePoly.gcdWith_eq.{u} {F : Type u} [DecidableEq F]
  [Lean.Grind.Field F] (plan : Hex.DensePoly.MulPlan F)
  (p q : Hex.DensePoly F) : Hex.DensePoly.gcdWith plan p q = p.gcd q
Hex.DensePoly.gcdWith_eq.{u} {F : Type u}
  [DecidableEq F] [Lean.Grind.Field F]
  (plan : Hex.DensePoly.MulPlan F)
  (p q : Hex.DensePoly F) :
  Hex.DensePoly.gcdWith plan p q = p.gcd q

Hex.DensePoly.xgcdLeftWith returns only the gcd and the coefficient of a, as HexPoly's xgcdLeft does.