Foundation Dalembert Inevitability P Symmetric From F Symmetric

/-- If the combiner `P` is symmetric and `F` is multiplicatively consistent with `P`,
then `F` is reciprocal-symmetric: `F(x) = F(1/x)` for all `x>0`. -/
theorem F_symmetric_of_P_symmetric (F : ℝ → ℝ) (P : ℝ → ℝ → ℝ)
    (hCons : HasMultiplicativeConsistency F P)
    (hSymP : ∀ u v, P u v = P v u) :
    IsSymmetric F := by
  intro x hx
  have h := F_div_swap_of_P_symmetric F P hCons hSymP x 1 hx one_pos
  simpa [div_one] using h

/-- **THEOREM: The consistency requirement forces the unique bilinear family.**

Given:
1. F : ℝ₊ → ℝ is a cost functional
2. F is symmetric: F(x) = F(1/x)
3. F is normalized: F(1) = 0
4. F has multiplicative consistency: F(xy) + F(x/y) = P(F(x), F(y)) for some **symmetric quadratic polynomial** P
5. F is non-trivial (not constant 0)

Then:
P(u, v) = 2u + 2v + c*u*v for some constant c.

This means F satisfies the generalized d'Alembert equation.
If we choose the canonical cost normalization c = 2, we recover the RCL. -/
theorem bilinear_family_forced (F : ℝ → ℝ) (P : ℝ → ℝ → ℝ)
    (hNorm : IsNormalized F)
    (hCons : HasMultiplicativeConsistency F P)
    (hPoly : ∃ (a b c d e f : ℝ), ∀ u v, P u v = a + b*u + c*v + d*u*v + e*u^2 + f*v^2)
    (hSymP : ∀ u v, P u v = P v u) -- Explicit symmetry of P
    (hNonTriv : ∃ x : ℝ, 0 < x ∧ F x ≠ 0)
    (hCont : ContinuousOn F (Set.Ioi 0)) -- Regularity: F is continuous on (0, ∞)
    : ∃ c : ℝ, (∀ u v, P u v = 2*u + 2*v + c*u*v) ∧
               (c = 2 → ∀ u v, P u v = 2*u + 2*v + 2*u*v) := by
  -- Derived reciprocity from symmetry of P
  have hSym : IsSymmetric F := F_symmetric_of_P_symmetric F P hCons hSymP
  -- Step 1: Normalization forces P(0, v) = 2v
  have hP0 : ∀ y : ℝ, 0 < y → P 0 (F y) = 2 * F y :=
    symmetry_and_normalization_constrain_P F P hSym hNorm hCons

  -- Use the polynomial form lemma
  -- We need to satisfy the hypotheses of `polynomial_form_forced`.
  -- `hNorm0`: ∀ v, P 0 v = 2 * v.
  -- We only have `P 0 (F y) = 2 F y`.
  -- However, since P is a polynomial and F is non-trivial (has range with at least 0 and some non-zero value),
  -- we can determine the coefficients.
  -- P(0, v) = a + c*v + f*v^2.
  -- P(0, 0) = a = 2*0 = 0 (from F(1)=0).
  -- P(0, F y) = c*(F y) + f*(F y)^2 = 2*(F y).
  -- This holds for y=1 (0=0) and some y where F y ≠ 0.
  -- If we only have two points, we can't uniquely determine a quadratic.
  -- But wait, `polynomial_form_forced` derived `a=0, c=2, f=0`.
  -- Let's reproduce that logic but being careful about the domain.

  obtain ⟨a, b, c, d, e, f, hP⟩ := hPoly

  -- 1. a = 0
  have ha : a = 0 := by
    have hCons1 := hCons 1 1 one_pos one_pos
    simp only [one_mul, one_div] at hCons1
    -- hCons1 : F 1 + F 1⁻¹ = P (F 1) (F 1)
    -- inv_one : 1⁻¹ = 1
    rw [inv_one, hNorm] at hCons1
    -- hCons1 : 0 + 0 = P 0 0
    simp only [add_zero] at hCons1
    -- hCons1 : 0 = P 0 0
    rw [hP 0 0] at hCons1
    simp at hCons1
    exact hCons1.symm

  -- 2. From hSymP: P(u,v) = P(v,u)
  -- a + bu + cv + duv + eu^2 + fv^2 = a + bv + cu + duv + ev^2 + fu^2
  -- (b-c)u + (c-b)v + (e-f)u^2 + (f-e)v^2 = 0
  -- This implies b=c and e=f.
  have hb_c : b = c := by
    have h1 := hSymP 1 0
    rw [hP 1 0, hP 0 1] at h1
    -- h1 : a + b*1 + c*0 + d*0 + e*1 + f*0 = a + b*0 + c*1 + d*0 + e*0 + f*1
    -- i.e., a + b + e = a + c + f
    -- Using ha: a = 0, we get b + e = c + f
    simp only [ha, mul_zero, mul_one, add_zero, zero_add] at h1
    -- We need another equation to separate b, e, c, f
    have h2 := hSymP 2 0
    rw [hP 2 0, hP 0 2] at h2
    simp only [ha, mul_zero, add_zero, zero_add] at h2
    -- h1: b + e = c + f
    -- h2: 2b + 4e = 2c + 4f
    -- From h2: b + 2e = c + 2f
    -- Subtracting h1: e = f
    -- So b = c
    linarith
  have he_f : e = f := by
    have h1 := hSymP 1 0
    have h2 := hSymP 2 0
    rw [hP 1 0, hP 0 1] at h1
    rw [hP 2 0, hP 0 2] at h2
    simp only [ha, mul_zero, mul_one, add_zero, zero_add] at h1 h2
    linarith

  -- Now P(0, v) = c*v + f*v^2 (using a=0, b=c, e=f).
  -- And P(0, F y) = 2 * F y.
  -- So c*(F y) + f*(F y)^2 = 2*(F y).
  -- (c - 2)*(F y) + f*(F y)^2 = 0.
  -- This must hold for all y > 0.
  -- Since F is non-trivial, there exists y such that F y ≠ 0.
  obtain ⟨y0, hy0_pos, hy0_ne⟩ := hNonTriv
  have hc_2 : c = 2 ∧ f = 0 := by
    -- Let k = F y0 (a nonzero value in the range).
    let k : ℝ := F y0
    have hk_ne : k ≠ 0 := by
      -- hy0_ne : F y0 ≠ 0
      simpa [k] using hy0_ne

    -- The polynomial identity on the range: (c - 2) * F(y) + f * (F(y))^2 = 0.
    have poly_identity : ∀ y : ℝ, 0 < y → (c - 2) * (F y) + f * (F y)^2 = 0 := by
      intro y hy
      have h := hP0 y hy
      rw [hP 0 (F y)] at h
      simp [ha, hb_c, he_f] at h
      linarith

    -- Use IVT to find y1 with F y1 = k/2 (since F(1)=0 and F(y0)=k).
    have hF1 : F 1 = 0 := hNorm
    have hInterval_pos : Set.uIcc 1 y0 ⊆ Set.Ioi 0 := by
      intro x hx
      rcases hx with ⟨hx_lo, _hx_hi⟩
      have hmin_pos : 0 < min 1 y0 := lt_min one_pos hy0_pos
      exact lt_of_lt_of_le hmin_pos hx_lo
    have hContInterval : ContinuousOn F (Set.uIcc 1 y0) :=
      hCont.mono hInterval_pos
    have h1_mem : 1 ∈ Set.uIcc 1 y0 := Set.left_mem_uIcc
    have hy0_mem : y0 ∈ Set.uIcc 1 y0 := Set.right_mem_uIcc
    have hk2_in_image : k / 2 ∈ F '' Set.uIcc 1 y0 := by
      have hPreconn := isPreconnected_uIcc (a := 1) (b := y0)
      by_cases hk : 0 ≤ k
      · -- monotone direction: 0 ≤ k, so k/2 ∈ Icc 0 k
        have hIVT : Set.Icc 0 k ⊆ F '' Set.uIcc 1 y0 := by
          simpa [hF1, k] using hPreconn.intermediate_value h1_mem hy0_mem hContInterval
        have hk2_between : k / 2 ∈ Set.Icc 0 k := by
          constructor <;> linarith
        exact hIVT hk2_between
      · -- reverse direction: k < 0, so k/2 ∈ Icc k 0
        have hk_lt : k < 0 := lt_of_not_ge hk
        have hIVT : Set.Icc k 0 ⊆ F '' Set.uIcc 1 y0 := by
          simpa [hF1, k] using hPreconn.intermediate_value hy0_mem h1_mem hContInterval
        have hk2_between : k / 2 ∈ Set.Icc k 0 := by
          constructor <;> linarith
        exact hIVT hk2_between
    obtain ⟨y1, hy1_mem, hFy1⟩ := hk2_in_image
    have hy1_pos : 0 < y1 := hInterval_pos hy1_mem

    -- Evaluate the polynomial identity at y0 and y1, then solve for c and f.
    have h_y0 : (c - 2) * k + f * k^2 = 0 := by
      have h := poly_identity y0 hy0_pos
      simpa [k] using h
    have h_y1 : (c - 2) * (k / 2) + f * (k / 2)^2 = 0 := by
      have h := poly_identity y1 hy1_pos
      -- rewrite F y1 = k/2
      simpa [hFy1, k] using h

    -- Multiply the y1 equation by 4 to align it with the y0 equation.
    have h_y1_4 : 2 * (c - 2) * k + f * k^2 = 0 := by
      have h' := congrArg (fun z => 4 * z) h_y1
      -- simplify 4*(...) and 4*0
      ring_nf at h'
      -- `ring_nf` chooses its own normal form; bridge to our preferred one.
      have hrew : c * k * 2 - k * 4 + k ^ 2 * f = 2 * (c - 2) * k + f * k ^ 2 := by ring
      -- h' : c*k*2 - k*4 + k^2*f = 0
      calc
        2 * (c - 2) * k + f * k ^ 2
            = c * k * 2 - k * 4 + k ^ 2 * f := by simpa [hrew] using (Eq.symm hrew)
        _ = 0 := h'

    -- Subtract to get (c - 2) * k = 0, hence c = 2 (since k ≠ 0).
    have hk_mul : (c - 2) * k = 0 := by
      linarith [h_y0, h_y1_4]
    have hc : c = 2 := by
      rcases mul_eq_zero.mp hk_mul with hc0 | hk0
      · linarith
      · exact False.elim (hk_ne hk0)

    -- Plug back to get f = 0.
    have hf : f = 0 := by
      have hk2_ne : k^2 ≠ 0 := pow_ne_zero 2 hk_ne
      have hfk2 : f * k^2 = 0 := by
        -- from h_y0 with c=2
        simpa [hc] using h_y0
      rcases mul_eq_zero.mp hfk2 with hf0 | hk20
      · exact hf0
      · exact False.elim (hk2_ne hk20)

    exact ⟨hc, hf⟩

  have hc : c = 2 := hc_2.1
  have hf : f = 0 := hc_2.2
  have hb : b = 2 := by rw [hb_c, hc]
  have he : e = 0 := by rw [he_f, hf]

  -- So P(u, v) = 2u + 2v + d*u*v.
  use d
  constructor
  · intro u v
    rw [hP, ha, hb, hc, he, hf]
    ring
  · intro hd u v
    rw [hP, ha, hb, hc, he, hf, hd]
    ring