8.3. Division by Newton iteration🔗

Division with remainder reduces to multiplication. Write rev f for the polynomial whose coefficients are those of f in reverse order. If a = q b + r with a of degree m, b of degree n and r of degree less than n, then reversing gives rev a ≡ rev q · rev b modulo xᵐ⁻ⁿ⁺¹. The constant term of rev b is the leading coefficient of b. When it is invertible, rev b has an inverse as a power series, and

rev q = rev a · (rev b)⁻¹ mod xᵐ⁻ⁿ⁺¹.

The remainder is then a - q b. Newton's method computes the inverse of a power series g: if h agrees with g⁻¹ to k terms, then h (2 - g h) agrees with it to 2k terms. Each step costs a few multiplications at the new precision, and the precisions double, so the whole inverse costs a constant times one multiplication at the final precision. Dividing a polynomial with 2n coefficients by one with n therefore costs a constant times one product of polynomials with n coefficients. Hex.DensePoly.reverseSeries reverses a polynomial into a series, and Hex.DensePoly.reciprocalWith runs Newton's method. It takes the inverse of the constant term as an argument, and the theorem Hex.DensePoly.reciprocalWith_eq says that it returns the same series as TSeries.invOfUnit from HexTruncatedSeries, which computes the inverse one coefficient at a time.

When many polynomials are divided by the same divisor, the inverse needs to be computed only once. A Hex.DensePoly.DivPlan stores a divisor and the inverse of its reversal to a fixed number of terms, its capacity, which bounds the length of the quotients it can produce. Hex.DensePoly.DivPlan.ofMonic builds one for a monic divisor over a commutative ring, and Hex.DensePoly.DivPlan.ofNonzero for a nonzero divisor over a field. Hex.DensePoly.DivPlan.divMod and Hex.DensePoly.DivPlan.mod take a proof that the quotient fits.

Modulo x² - x - 1, the power xᵏ reduces to Fₖ x + Fₖ₋₁, where Fₖ is the kth Fibonacci number:

def golden : DensePoly Int := #p[-1, -1, 1]

def goldenPlan : DivPlan Int :=
  DivPlan.ofMonic (karatsubaPlan 32) golden
    (by⊢ golden.Monic unfold Monic⊢ golden.leadingCoeff = 1; decide +kernelAll goals completed! 🐙) (by⊢ golden ≠ 0 decideAll goals completed! 🐙) 128

def reduce (p : DensePoly Int) :
    Option (DensePoly Int) :=
  if h : quotientLength p goldenPlan.divisor ≤
      goldenPlan.capacity then
    some (goldenPlan.mod p h)
  else none

some #p[218922995834555169026, 354224848179261915075]#eval reduce (monomial 100 1)

some #p[218922995834555169026, 354224848179261915075]

For a single division over a field, Hex.DensePoly.divModWith builds the plan and uses it once. It returns exactly the quotient and remainder of HexPoly's long division, including its conventions for a zero divisor and for a divisor longer than the dividend:

theorem

Hex.DensePoly.divModWith_eq.{u} {F : Type u} [DecidableEq F]
  [Lean.Grind.Field F] (mul : Hex.DensePoly.MulPlan F)
  (p q : Hex.DensePoly F) :
  Hex.DensePoly.divModWith mul p q = p.divMod q
Hex.DensePoly.divModWith_eq.{u}
  {F : Type u} [DecidableEq F]
  [Lean.Grind.Field F]
  (mul : Hex.DensePoly.MulPlan F)
  (p q : Hex.DensePoly F) :
  Hex.DensePoly.divModWith mul p q =
    p.divMod q

Planned division has the exact reference quotient and remainder.

When the divisor or the quotient has at most eight coefficients, long division is faster, and compiled code for divModWith uses it. In fact divModWith is defined to be HexPoly's divMod, and a @[csimp] lemma, which proves the two functions equal, tells the compiler to run the Newton version instead. So a proof that evaluates divModWith in the kernel, for example by decide, sees long division. The gcd functions and multipoint evaluation below are set up in the same way.