8.5. Many points at once🔗

Evaluating a polynomial of degree less than n at n points by Horner's rule takes n² operations. The remainder of f modulo x - a is f(a), and the remainders of f modulo many polynomials can be computed together with a product tree. Its leaves are the moduli m₁, …, mₙ, and each node holds the product of the two nodes below it, so the root is m₁ ⋯ mₙ. The remainder tree reduces f modulo the root, then reduces that remainder modulo the two children of the root, and so on down to the leaves. With fast multiplication and division, for n linear leaves both trees take O(M(n) log n) operations.

Hex.DensePoly.ProductTree.build builds a product tree from an array of polynomials. It keeps only the leaves and the root, and recomputes the other levels on request with Hex.DensePoly.ProductTree.level?. Hex.DensePoly.RemainderTree.build builds a remainder tree from monic leaves, Hex.DensePoly.MonicLeaf, and stores a division plan at every node. Hex.DensePoly.RemainderTree.remainders? returns the remainders of a polynomial p modulo every leaf, in order; it succeeds whenever p.size ≤ tree.rootDegree + tree.capacity, where rootDegree is the sum of the degrees of the leaves and capacity is chosen when the tree is built.

Hex.DensePoly.EvalPlan.build builds the remainder tree for the leaves x - aᵢ at a list of points, with capacity for polynomials with as many coefficients as there are points. Hex.DensePoly.EvalPlan.eval then evaluates a polynomial at all the points; a longer polynomial is evaluated at each point by Horner's rule.

Hex.DensePoly.EvalPlan.get_eval.{u} {R : Type u} [DecidableEq R]
  [Lean.Grind.CommRing R] (plan : Hex.DensePoly.EvalPlan R)
  (f : Hex.DensePoly R) (i : ℕ) (hi : i < plan.size) :
  (plan.eval f)[i] = f.eval plan.points[i]
Hex.DensePoly.EvalPlan.get_eval.{u}
  {R : Type u} [DecidableEq R]
  [Lean.Grind.CommRing R]
  (plan : Hex.DensePoly.EvalPlan R)
  (f : Hex.DensePoly R) (i : ℕ)
  (hi : i < plan.size) :
  (plan.eval f)[i] = f.eval plan.points[i]

Interpolation inverts evaluation, and needs a field. Given distinct points a₁, …, aₙ and values y₁, …, yₙ, let m = (x - a₁) ⋯ (x - aₙ). Lagrange's formula for the unique polynomial of degree less than n with f(aᵢ) = yᵢ is

f = Σᵢ yᵢ / m'(aᵢ) · m / (x - aᵢ).

Hex.DensePoly.InterpPlan.build? builds the product tree for m, evaluates the derivative m' at all the points with the remainder tree, and stores the inverses 1 / m'(aᵢ). It returns none exactly when two points are equal. Hex.DensePoly.InterpPlan.interpolate? then sums Lagrange's formula up the tree: a node whose children have products m_L and m_R, and partial sums f_L and f_R, gets f_L m_R + f_R m_L. It returns none exactly when the number of values differs from the number of points.

The sums 0, 1, 5, 14 of the squares up to 0², 1², 2², 3² determine the cubic polynomial n(n + 1)(2n + 1) / 6 that gives the sum of the squares up to n²:

some #p[0, (1 : Rat)/6, (1 : Rat)/2, (1 : Rat)/3]

Evaluating at 4, 5, 6, 7 continues the sequence:

some #[30, 55, 91, 140]

Hex.DensePoly.InterpPlan.interpolate?_sound.{u} {F : Type u}
  [DecidableEq F] [Lean.Grind.Field F]
  (plan : Hex.DensePoly.InterpPlan F) (values : Array F)
  (p : Hex.DensePoly F) (hsome : plan.interpolate? values = some p) :
  p.size ≤ plan.size ∧
    ∀ (i : ℕ) (hi : i < plan.size),
      p.eval plan.points[i] = values.getD i 0
Hex.DensePoly.InterpPlan.interpolate?_sound.{u}
  {F : Type u} [DecidableEq F]
  [Lean.Grind.Field F]
  (plan : Hex.DensePoly.InterpPlan F)
  (values : Array F) (p : Hex.DensePoly F)
  (hsome :
    plan.interpolate? values = some p) :
  p.size ≤ plan.size ∧
    ∀ (i : ℕ) (hi : i < plan.size),
      p.eval plan.points[i] =
        values.getD i 0

Hex.DensePoly.InterpPlan.interpolate?_unique.{u} {F : Type u}
  [DecidableEq F] [Lean.Grind.Field F]
  (plan : Hex.DensePoly.InterpPlan F) (values : Array F)
  (p q : Hex.DensePoly F) (hsome : plan.interpolate? values = some p)
  (hqsize : q.size ≤ plan.size)
  (hqvalues :
    ∀ (i : ℕ) (hi : i < plan.size),
      q.eval plan.points[i] = values.getD i 0) :
  q = p
Hex.DensePoly.InterpPlan.interpolate?_unique.{u}
  {F : Type u} [DecidableEq F]
  [Lean.Grind.Field F]
  (plan : Hex.DensePoly.InterpPlan F)
  (values : Array F)
  (p q : Hex.DensePoly F)
  (hsome :
    plan.interpolate? values = some p)
  (hqsize : q.size ≤ plan.size)
  (hqvalues :
    ∀ (i : ℕ) (hi : i < plan.size),
      q.eval plan.points[i] =
        values.getD i 0) :
  q = p