Every scalar real polynomial minimal graph in dimension eight is affine

/- Every scalar real polynomial minimal graph in dimension eight is affine.
   Submission body for Math30Catalog.source29.
   The challenge supplies its imports and the surrounding Bounty namespace. -/

/- Component: LeadingReduction.lean -/

/-!
# All-degree leading-component reduction for polynomial minimal graphs

Extracted and adapted from the PolynomialBasics component of
`QuarticMinimalGraphs.lean` (the Math30Catalog.source28 submission body).
The statements below make no upper bound assumption on polynomial degree.
This is a proved dependency toward Math30Catalog.source29, not a proof of that target.
-/

/- Component: PolynomialBasics -/
section ComponentPolynomialBasics

noncomputable section
open scoped BigOperators
open MvPolynomial

namespace QuarticPolynomialBasics

variable {σ : Type*}


lemma totalDegree_pderiv_le_of_totalDegree_le (P : MvPolynomial σ ℝ)
    (i : σ) (d : ℕ) (h : P.totalDegree ≤ d + 1) :
    (pderiv i P).totalDegree ≤ d := by
  classical
  apply Finset.sup_le
  intro m hm
  change m.degree ≤ d
  have hc : P.coeff (m + Finsupp.single i 1) ≠ 0 := by
    have hh := mem_support_iff.mp hm
    rw [coeff_pderiv] at hh
    exact (mul_ne_zero_iff.mp hh).1
  have hd := le_totalDegree (mem_support_iff.mpr hc)
  change (m + Finsupp.single i 1).degree ≤ P.totalDegree at hd
  simp only [map_add, Finsupp.degree_single] at hd
  omega

lemma totalDegree_pderiv_le (P : MvPolynomial σ ℝ) (i : σ) :
    (pderiv i P).totalDegree ≤ P.totalDegree - 1 :=
  totalDegree_pderiv_le_of_totalDegree_le P i _ (by omega)

lemma pderiv_homogeneousComponent (P : MvPolynomial σ ℝ) (i : σ) (d : ℕ) :
    pderiv i (homogeneousComponent (d + 1) P) =
      homogeneousComponent d (pderiv i P) := by
  classical
  ext m
  simp only [coeff_pderiv, coeff_homogeneousComponent, map_add, Finsupp.degree_single]
  by_cases h : m.degree = d <;> simp [h]


lemma totalDegree_sub_homogeneousComponent_le (P : MvPolynomial σ ℝ)
    (d : ℕ) (h : P.totalDegree ≤ d) :
    (P - homogeneousComponent d P).totalDegree ≤ d - 1 := by
  classical
  apply Finset.sup_le
  intro m hm
  change m.degree ≤ d - 1
  have hc := mem_support_iff.mp hm
  simp only [coeff_sub, coeff_homogeneousComponent] at hc
  split_ifs at hc with hd
  · exact False.elim (hc (sub_self _))
  · simp only [sub_zero] at hc
    have hh := le_totalDegree (mem_support_iff.mpr hc)
    change m.degree ≤ P.totalDegree at hh
    omega

lemma top_homogeneousComponent_ne_zero (P : MvPolynomial σ ℝ)
    (h : 0 < P.totalDegree) : homogeneousComponent P.totalDegree P ≠ 0 := by
  intro hh
  have hdeg := totalDegree_sub_homogeneousComponent_le P P.totalDegree le_rfl
  rw [hh, sub_zero] at hdeg
  omega


lemma homogeneousComponent_mul_top (P Q : MvPolynomial σ ℝ) (d e : ℕ)
    (hP : P.totalDegree ≤ d) (hQ : Q.totalDegree ≤ e) :
    homogeneousComponent (d + e) (P * Q) =
      homogeneousComponent d P * homogeneousComponent e Q := by
  by_cases hd : d = 0
  · subst d
    have hp0 : P.totalDegree = 0 := Nat.eq_zero_of_le_zero hP
    rw [(totalDegree_eq_zero_iff_eq_C.mp hp0)]
    simp [homogeneousComponent_C_mul, homogeneousComponent_zero]
  by_cases he : e = 0
  · subst e
    have hq0 : Q.totalDegree = 0 := Nat.eq_zero_of_le_zero hQ
    rw [(totalDegree_eq_zero_iff_eq_C.mp hq0)]
    rw [mul_comm P (C (Q.coeff 0)), homogeneousComponent_C_mul]
    simp [homogeneousComponent_zero, mul_comm]
  have hdp := totalDegree_sub_homogeneousComponent_le P d hP
  have hdq := totalDegree_sub_homogeneousComponent_le Q e hQ
  have htp : (homogeneousComponent d P).totalDegree ≤ d :=
    (homogeneousComponent_isHomogeneous d P).totalDegree_le
  have hz1 : homogeneousComponent (d + e) ((P - homogeneousComponent d P) * Q) = 0 := by
    apply homogeneousComponent_eq_zero
    have hh := totalDegree_mul (P - homogeneousComponent d P) Q
    omega
  have hz2 : homogeneousComponent (d + e)
      (homogeneousComponent d P * (Q - homogeneousComponent e Q)) = 0 := by
    apply homogeneousComponent_eq_zero
    have hh := totalDegree_mul (homogeneousComponent d P) (Q - homogeneousComponent e Q)
    omega
  have heq : P * Q = homogeneousComponent d P * homogeneousComponent e Q +
      ((P - homogeneousComponent d P) * Q +
        homogeneousComponent d P * (Q - homogeneousComponent e Q)) := by ring
  rw [heq, map_add, map_add, hz1, hz2, add_zero, add_zero]
  exact homogeneousComponent_eq_self
    ((homogeneousComponent_isHomogeneous d P).mul (homogeneousComponent_isHomogeneous e Q))

/-- The polynomial numerator of the scalar minimal graph equation. This is
definitionally equal to the corresponding operator in the challenge. -/
def minimalGraphOperator {n : ℕ} (P : MvPolynomial (Fin n) ℝ) :
    MvPolynomial (Fin n) ℝ :=
  (1 + ∑ i : Fin n, (pderiv i P) ^ 2) *
      (∑ i : Fin n, pderiv i (pderiv i P)) -
    ∑ i : Fin n, ∑ j : Fin n,
      pderiv i P * pderiv j P * pderiv i (pderiv j P)


def polynomialLaplacian {n : ℕ} (P : MvPolynomial (Fin n) ℝ) :
    MvPolynomial (Fin n) ℝ := ∑ i : Fin n, pderiv i (pderiv i P)

def cubicOperator {n : ℕ} (P : MvPolynomial (Fin n) ℝ) :
    MvPolynomial (Fin n) ℝ :=
  (∑ i : Fin n, (pderiv i P) ^ 2) * polynomialLaplacian P -
    ∑ i : Fin n, ∑ j : Fin n,
      pderiv i P * pderiv j P * pderiv i (pderiv j P)

lemma minimalGraphOperator_eq_laplacian_add_cubic {n : ℕ}
    (P : MvPolynomial (Fin n) ℝ) :
    minimalGraphOperator P = polynomialLaplacian P + cubicOperator P := by
  unfold minimalGraphOperator cubicOperator polynomialLaplacian
  ring

lemma totalDegree_second_partial_le (P : MvPolynomial σ ℝ) (d : ℕ)
    (h : P.totalDegree ≤ d + 2) (i j : σ) :
    (pderiv i (pderiv j P)).totalDegree ≤ d := by
  apply totalDegree_pderiv_le_of_totalDegree_le
  exact totalDegree_pderiv_le_of_totalDegree_le P j (d + 1) (by omega)

lemma homogeneousComponent_second_partial (P : MvPolynomial σ ℝ) (d : ℕ)
    (i j : σ) : homogeneousComponent d (pderiv i (pderiv j P)) =
      pderiv i (pderiv j (homogeneousComponent (d + 2) P)) := by
  rw [← pderiv_homogeneousComponent, ← pderiv_homogeneousComponent]

lemma homogeneousComponent_derivative_triple (P : MvPolynomial σ ℝ) (d : ℕ)
    (h : P.totalDegree ≤ d + 2) (i j k l : σ) :
    homogeneousComponent (3 * d + 2)
      (pderiv i P * pderiv j P * pderiv k (pderiv l P)) =
      pderiv i (homogeneousComponent (d + 2) P) *
        pderiv j (homogeneousComponent (d + 2) P) *
        pderiv k (pderiv l (homogeneousComponent (d + 2) P)) := by
  have hgi : (pderiv i P).totalDegree ≤ d + 1 :=
    totalDegree_pderiv_le_of_totalDegree_le P i (d + 1) (by omega)
  have hgj : (pderiv j P).totalDegree ≤ d + 1 :=
    totalDegree_pderiv_le_of_totalDegree_le P j (d + 1) (by omega)
  have hpair : (pderiv i P * pderiv j P).totalDegree ≤ (d + 1) + (d + 1) :=
    (totalDegree_mul _ _).trans (Nat.add_le_add hgi hgj)
  have hkl := totalDegree_second_partial_le P d h k l
  rw [show 3 * d + 2 = ((d + 1) + (d + 1)) + d by omega,
    homogeneousComponent_mul_top _ _ _ _ hpair hkl,
    homogeneousComponent_mul_top _ _ _ _ hgi hgj,
    ← pderiv_homogeneousComponent, ← pderiv_homogeneousComponent,
    homogeneousComponent_second_partial]

lemma homogeneousComponent_cubicOperator_top {n : ℕ}
    (P : MvPolynomial (Fin n) ℝ) (d : ℕ) (h : P.totalDegree ≤ d + 2) :
    homogeneousComponent (3 * d + 2) (cubicOperator P) =
      cubicOperator (homogeneousComponent (d + 2) P) := by
  unfold cubicOperator polynomialLaplacian
  simp only [Finset.sum_mul, Finset.mul_sum, map_sub, map_sum, pow_two]
  simp_rw [homogeneousComponent_derivative_triple P d h]

lemma homogeneousComponent_minimalGraphOperator_top {n : ℕ}
    (P : MvPolynomial (Fin n) ℝ) (d : ℕ) (h : P.totalDegree ≤ d + 2) :
    homogeneousComponent (3 * d + 2) (minimalGraphOperator P) =
      cubicOperator (homogeneousComponent (d + 2) P) := by
  rw [minimalGraphOperator_eq_laplacian_add_cubic, map_add,
    homogeneousComponent_cubicOperator_top P d h]
  have hlap : (polynomialLaplacian P).totalDegree ≤ d := by
    apply totalDegree_finsetSum_le
    intro i hi
    exact totalDegree_second_partial_le P d h i i
  have hz : homogeneousComponent (3 * d + 2) (polynomialLaplacian P) = 0 := by
    apply homogeneousComponent_eq_zero
    omega
  rw [hz, zero_add]

lemma cubicOperator_top_eq_zero {n : ℕ} (P : MvPolynomial (Fin n) ℝ)
    (m : ℕ) (hm : 2 ≤ m) (hdeg : P.totalDegree ≤ m)
    (h : minimalGraphOperator P = 0) :
    cubicOperator (homogeneousComponent m P) = 0 := by
  obtain ⟨d, rfl⟩ := Nat.exists_eq_add_of_le hm
  have ht := homogeneousComponent_minimalGraphOperator_top P d (by omega)
  rw [h, map_zero] at ht
  simpa [Nat.add_comm] using ht.symm

/-- Every nonaffine polynomial minimal graph has a nonzero leading homogeneous
component annihilated by the cubic part of the minimal graph operator. -/
theorem leading_component_data {n : ℕ} (P : MvPolynomial (Fin n) ℝ)
    (hmin : Math30.MinimalGraphs.SatisfiesMinimalGraphEquation P)
    (hdeg : 2 ≤ P.totalDegree) :
    let F := homogeneousComponent P.totalDegree P
    F ≠ 0 ∧ F.IsHomogeneous P.totalDegree ∧ cubicOperator F = 0 := by
  refine ⟨top_homogeneousComponent_ne_zero P (by omega),
    homogeneousComponent_isHomogeneous _ _, ?_⟩
  exact cubicOperator_top_eq_zero P P.totalDegree hdeg le_rfl hmin

/-- The strict degree interval appearing in dimension eight contains no integer. -/
theorem no_integer_degree_between_three_and_four (d : ℕ)
    (hlower : (3 : ℝ) < d) (hupper : (d : ℝ) < 4) : False := by
  have hlower' : 3 < d := by exact_mod_cast hlower
  have hupper' : d < 4 := by exact_mod_cast hupper
  omega

end QuarticPolynomialBasics
end
end ComponentPolynomialBasics

/- Component: HomogeneousFactor.lean -/

/-!
# Homogeneous factors of a nonzero homogeneous polynomial

The proof encodes total degree by a new univariate parameter. The largest and
smallest degrees of a nonzero product add, so a product supported in one degree
forces each factor to be supported in one degree as well.
-/

namespace Source29HomogeneousFactor

open MvPolynomial
open scoped BigOperators

variable {σ R : Type*} [CommRing R]

noncomputable def degreeEncoding :
    MvPolynomial σ R →+* Polynomial (MvPolynomial σ R) :=
  eval₂Hom (Polynomial.C.comp C)
    (fun i => Polynomial.C (X i) * Polynomial.X)

lemma degreeEncoding_monomial (e : σ →₀ ℕ) (c : R) :
    degreeEncoding (monomial e c) = Polynomial.monomial e.degree (monomial e c) := by
  classical
  rw [degreeEncoding, eval₂Hom_monomial]
  simp only [RingHom.comp_apply, Finsupp.prod, mul_pow,
    Finset.prod_mul_distrib, ← Polynomial.C_pow, ← map_prod,
    Finset.prod_pow_eq_pow_sum, ← Polynomial.C_mul_X_pow_eq_monomial]
  rw [← mul_assoc, ← Polynomial.C_mul]
  change Polynomial.C (C c * e.prod (fun i k => X i ^ k)) * Polynomial.X ^ e.degree =
    Polynomial.C (monomial e c) * Polynomial.X ^ e.degree
  rw [← monomial_eq]

lemma coeff_degreeEncoding (P : MvPolynomial σ R) (d : ℕ) :
    (degreeEncoding P).coeff d = homogeneousComponent d P := by
  classical
  induction P using MvPolynomial.induction_on' with
  | monomial e c =>
    rw [degreeEncoding_monomial, Polynomial.coeff_monomial,
      homogeneousComponent_of_mem (isHomogeneous_monomial c rfl)]
    simp [eq_comm]
  | add P Q hP hQ => simp only [map_add, Polynomial.coeff_add, hP, hQ]

lemma degreeEncoding_ne_zero {P : MvPolynomial σ R} (hP : P ≠ 0) :
    degreeEncoding P ≠ 0 := by
  intro hzero
  apply hP
  rw [← sum_homogeneousComponent P]
  apply Finset.sum_eq_zero
  intro d hd
  rw [← coeff_degreeEncoding, hzero, Polynomial.coeff_zero]

lemma degreeEncoding_of_homogeneous {P : MvPolynomial σ R} {m : ℕ}
    (hP : P.IsHomogeneous m) :
    degreeEncoding P = Polynomial.monomial m P := by
  ext d
  rw [coeff_degreeEncoding, homogeneousComponent_of_mem hP,
    Polynomial.coeff_monomial]
  simp [eq_comm]

lemma homogeneous_of_encoding_equal_degrees {P : MvPolynomial σ R}
    (hdegree : (degreeEncoding P).natDegree = (degreeEncoding P).natTrailingDegree) :
    P.IsHomogeneous (degreeEncoding P).natDegree := by
  classical
  intro e he
  change (Finsupp.weight fun _ : σ => (1 : ℕ)) e = _
  rw [← Finsupp.degree_eq_weight_one]
  have hcoeff : (degreeEncoding P).coeff e.degree ≠ 0 := by
    rw [coeff_degreeEncoding]
    intro hz
    have h := congrArg (fun Q : MvPolynomial σ R => Q.coeff e) hz
    apply he
    simpa [coeff_homogeneousComponent] using h
  have hlo := Polynomial.natTrailingDegree_le_of_ne_zero hcoeff
  have hhi := Polynomial.le_natDegree_of_ne_zero hcoeff
  omega

variable [IsDomain R]

/-- Every factor of a nonzero homogeneous product is homogeneous. -/
theorem left_factor_isHomogeneous (P Q : MvPolynomial σ R) {m : ℕ}
    (hne : P * Q ≠ 0) (hhom : (P * Q).IsHomogeneous m) :
    P.IsHomogeneous P.totalDegree := by
  have hP : P ≠ 0 := left_ne_zero_of_mul hne
  have hQ : Q ≠ 0 := right_ne_zero_of_mul hne
  have hA := degreeEncoding_ne_zero hP
  have hB := degreeEncoding_ne_zero hQ
  have heq : degreeEncoding P * degreeEncoding Q = Polynomial.monomial m (P * Q) := by
    rw [← map_mul]
    exact degreeEncoding_of_homogeneous hhom
  have htop := congrArg Polynomial.natDegree heq
  rw [Polynomial.natDegree_mul hA hB, Polynomial.natDegree_monomial_eq m hne] at htop
  have hbottom := congrArg Polynomial.natTrailingDegree heq
  rw [Polynomial.natTrailingDegree_mul hA hB,
    Polynomial.natTrailingDegree_monomial hne] at hbottom
  have hAle := Polynomial.natTrailingDegree_le_natDegree (degreeEncoding P)
  have hBle := Polynomial.natTrailingDegree_le_natDegree (degreeEncoding Q)
  have hsame : (degreeEncoding P).natDegree = (degreeEncoding P).natTrailingDegree := by omega
  have hhomP := homogeneous_of_encoding_equal_degrees hsame
  simpa only [hhomP.totalDegree hP] using hhomP

/-- A divisor of a nonzero homogeneous polynomial is homogeneous of its total degree. -/
theorem divisor_isHomogeneous (F h : MvPolynomial σ R) {m : ℕ}
    (hne : F ≠ 0) (hhom : F.IsHomogeneous m) (hdiv : h ∣ F) :
    h.IsHomogeneous h.totalDegree := by
  obtain ⟨g, rfl⟩ := hdiv
  exact left_factor_isHomogeneous h g hne hhom

end Source29HomogeneousFactor

/- Component: FactorInheritance.lean -/

namespace Source29FactorInheritance
open scoped BigOperators
open MvPolynomial

section Algebra
variable {R ι : Type*} [CommRing R] [Fintype ι]

noncomputable def form (g : ι → R) (U : ι → ι → R) : R :=
  ∑ i, ∑ j, (g i ^ 2 * U j j - g i * g j * U i j)

lemma form_scaled (s t : R) (g : ι → R) (U : ι → ι → R) :
    form (fun i => s * g i) (fun i j => t * U i j) = s ^ 2 * t * form g U := by
  simp only [form, Finset.mul_sum]
  apply Finset.sum_congr rfl
  intro i hi
  apply Finset.sum_congr rfl
  intro j hj
  ring

noncomputable def remainderPair (t k v : R) (g b : ι → R) (U V : ι → ι → R)
    (i j : ι) : R :=
  (k^2*v^2*g i^2*V j j + 2*k^2*v*g i*b i*(v*U j j+2*g j*b j)
    + k*(k-1)*v*b i^2*g j^2 - k^2*v^2*g i*g j*V i j
    - k^2*v*(g i*b j+g j*b i)*(v*U i j+g i*b j+g j*b i)
    - k*(k-1)*v*b i*b j*g i*g j) +
  t * (2*k*v*g i*b i*V j j + k*b i^2*(v*U j j+2*g j*b j)
    - k*v*(g i*b j+g j*b i)*V i j
    - k*b i*b j*(v*U i j+g i*b j+g j*b i)) +
  t^2*(b i^2*V j j-b i*b j*V i j)

lemma form_normalized (t k v : R) (g b : ι → R) (U V : ι → ι → R) :
    ∃ r : R,
      form (fun i => k*v*g i+t*b i)
        (fun i j => k*(k-1)*v*g i*g j +
          t*k*(v*U i j+g i*b j+g j*b i)+t^2*V i j) =
      t*(k^3*v^3*form g U+t*r) := by
  classical
  let r := ∑ i, ∑ j, remainderPair t k v g b U V i j
  refine ⟨r, ?_⟩
  have hcross : (∑ i, ∑ j, (g i^2*g j*b j-g i*g j^2*b i)) = 0 := by
    simp only [Finset.sum_sub_distrib, sub_eq_zero]
    rw [Finset.sum_comm]
    apply Finset.sum_congr rfl
    intro i hi
    apply Finset.sum_congr rfl
    intro j hj
    ring
  calc
    _ = ∑ i, ∑ j, t * (k^3*v^3*(g i^2*U j j-g i*g j*U i j) +
        k^2*v^2*(g i^2*g j*b j-g i*g j^2*b i) +
        t*remainderPair t k v g b U V i j) := by
      unfold form
      apply Finset.sum_congr rfl
      intro i hi
      apply Finset.sum_congr rfl
      intro j hj
      dsimp [remainderPair]
      ring
    _ = t*(k^3*v^3*form g U + k^2*v^2*
        (∑ i, ∑ j, (g i^2*g j*b j-g i*g j^2*b i)) + t*r) := by
      simp only [r, form, mul_add, Finset.sum_add_distrib, ← Finset.mul_sum]
    _ = _ := by rw [hcross]; ring

end Algebra

section Polynomial
variable {R ι : Type*} [CommRing R] [Fintype ι]

noncomputable def cubicOperator (f : MvPolynomial ι R) : MvPolynomial ι R :=
  (∑ i, (pderiv i f)^2)*(∑ j, pderiv j (pderiv j f)) -
    ∑ i, ∑ j, pderiv i f*pderiv j f*pderiv i (pderiv j f)

lemma cubicOperator_eq_form (f : MvPolynomial ι R) :
    cubicOperator f = form (fun i => pderiv i f) (fun i j => pderiv i (pderiv j f)) := by
  simp only [cubicOperator, form, Finset.sum_sub_distrib, Finset.sum_mul,
    Finset.mul_sum]
  rw [Finset.sum_comm]

lemma scaled_pderiv_pow_mul (h v : MvPolynomial ι R) (a : ℕ) (i : ι) :
    h * pderiv i (h^a*v) =
      h^a * ((a : MvPolynomial ι R)*v*pderiv i h + h*pderiv i v) := by
  cases a with
  | zero => simp
  | succ a =>
    rw [pderiv_mul, pderiv_pow, Nat.succ_sub_one]
    simp only [pow_succ, Nat.cast_add, Nat.cast_one]
    ring

lemma scaled_pderiv_pow (h : MvPolynomial ι R) (a : ℕ) (i : ι) :
    h * pderiv i (h^a) = h^a * (a : MvPolynomial ι R) * pderiv i h := by
  simpa only [mul_one, pderiv_one, mul_zero, add_zero, mul_assoc] using
    scaled_pderiv_pow_mul h 1 a i

lemma scaled_second_pderiv_pow_mul (h v : MvPolynomial ι R) (a : ℕ) (i j : ι) :
    h^2 * pderiv i (pderiv j (h^a*v)) = h^a *
      ((a : MvPolynomial ι R)*((a : MvPolynomial ι R)-1)*v*pderiv i h*pderiv j h +
      h*(a : MvPolynomial ι R)*(v*pderiv i (pderiv j h) +
        pderiv i h*pderiv j v + pderiv j h*pderiv i v) +
      h^2*pderiv i (pderiv j v)) := by
  let F := h^a*v
  have hg : h*pderiv j F = h^a*((a : MvPolynomial ι R)*v*pderiv j h+
      h*pderiv j v) := scaled_pderiv_pow_mul h v a j
  have hd := congrArg (pderiv i) hg
  simp only [pderiv_mul, map_add, Derivation.map_natCast, zero_mul, zero_add] at hd
  have hp := scaled_pderiv_pow h a i
  change h^2*pderiv i (pderiv j F) = _
  linear_combination h * hd - (pderiv i h) * hg +
    ((a : MvPolynomial ι R)*v*pderiv j h + h*pderiv j v) * hp


/-- A denominator-free form of factor inheritance, valid for every multiplicity. -/
lemma scaled_cubicOperator_pow_mul (h v : MvPolynomial ι R) (a : ℕ) :
    ∃ r : MvPolynomial ι R,
      h^4 * cubicOperator (h^a*v) =
      (h^a)^3*h*((a : MvPolynomial ι R)^3*v^3*cubicOperator h+h*r) := by
  obtain ⟨r, hr⟩ := form_normalized h (a : MvPolynomial ι R) v
    (fun i => pderiv i h) (fun i => pderiv i v)
    (fun i j => pderiv i (pderiv j h)) (fun i j => pderiv i (pderiv j v))
  refine ⟨r, ?_⟩
  calc
    _ = form (fun i => h*pderiv i (h^a*v))
        (fun i j => h^2*pderiv i (pderiv j (h^a*v))) := by
      rw [form_scaled, ← cubicOperator_eq_form]
      ring
    _ = form
        (fun i => h^a*((a : MvPolynomial ι R)*v*pderiv i h+h*pderiv i v))
        (fun i j => h^a*((a : MvPolynomial ι R)*((a : MvPolynomial ι R)-1)*v*
          pderiv i h*pderiv j h + h*(a : MvPolynomial ι R)*
          (v*pderiv i (pderiv j h)+pderiv i h*pderiv j v+pderiv j h*pderiv i v)+
          h^2*pderiv i (pderiv j v))) := by
      simp_rw [scaled_pderiv_pow_mul, scaled_second_pderiv_pow_mul]
    _ = _ := by
      rw [form_scaled, hr, ← cubicOperator_eq_form]
      ring


/-- The exact multiplicity identity stated in the proof candidate. -/
lemma cubicOperator_pow_mul [IsDomain R] (h v : MvPolynomial ι R) (a : ℕ)
    (ha : 1 ≤ a) (hh : h ≠ 0) :
    ∃ r : MvPolynomial ι R,
      cubicOperator (h^a*v) =
      h^(3*a-3)*((a : MvPolynomial ι R)^3*v^3*cubicOperator h+h*r) := by
  obtain ⟨r, hr⟩ := scaled_cubicOperator_pow_mul h v a
  refine ⟨r, ?_⟩
  apply mul_left_cancel₀ (pow_ne_zero 4 hh)
  rw [hr]
  have hp : (h^a)^3*h = h^(3*a+1) := by
    rw [← pow_mul, ← pow_succ]
    congr 1
    omega
  rw [hp, ← mul_assoc, ← pow_add]
  congr 2
  omega

end Polynomial

/-- Every irreducible factor of a polynomial annihilated by the cubic operator
inherits the operator equation modulo that factor, with arbitrary multiplicity. -/
theorem factor_inheritance {ι : Type*} [Fintype ι]
    (h v : MvPolynomial ι ℂ) (a : ℕ) (ha : 1 ≤ a)
    (hh : Irreducible h) (hv : ¬h ∣ v)
    (hF : cubicOperator (h^a*v) = 0) : h ∣ cubicOperator h := by