/-
EEML-kernstelling, machinale verificatie (13-8-2026).

Stelling: laat P : ℝ → ℝ op de niet-negatieve reële as voldoen aan
  (i)  compositie:  P (s+t) = P s * P t          [axioma A2]
  (ii) codomein:    0 < P s ≤ 1                   [codomein-eis]
Dan is er een λ ≥ 0 met P s = exp (−λ s) voor alle s ≥ 0.

De monotonie wordt niet aangenomen: zij volgt uit het codomein
(P (s+t) ≤ P s omdat P t ≤ 1), dus er is geen regulariteits-
aanname nodig naast de twee axioma's.
-/
import Mathlib

open Real

namespace EEML

/-- n-voudige additiviteit op ℝ≥0. -/
lemma iterate_add (g : ℝ → ℝ)
    (hadd : ∀ s t : ℝ, 0 ≤ s → 0 ≤ t → g (s + t) = g s + g t)
    (s : ℝ) (hs : 0 ≤ s) :
    ∀ n : ℕ, g (n * s) = n * g s := by
  intro n
  induction n with
  | zero =>
      have h : g (0 + 0) = g 0 + g 0 := hadd 0 0 le_rfl le_rfl
      simp at h
      simpa using h
  | succ k ih =>
      have hk : (0:ℝ) ≤ (k:ℝ) * s := by positivity
      have harg : (((k+1 : ℕ)):ℝ) * s = (k:ℝ) * s + s := by push_cast; ring
      rw [harg, hadd _ _ hk hs, ih]
      push_cast; ring

/-- Additief + niet-negatief op ℝ≥0 ⇒ lineair op ℝ≥0
    (Ostrowski-variant: monotonie volgt uit de tekens). -/
lemma additive_nonneg_linear (g : ℝ → ℝ)
    (hadd : ∀ s t : ℝ, 0 ≤ s → 0 ≤ t → g (s + t) = g s + g t)
    (hnn : ∀ s : ℝ, 0 ≤ s → 0 ≤ g s) :
    ∀ s : ℝ, 0 ≤ s → g s = g 1 * s := by
  -- monotonie
  have hmono : ∀ a b : ℝ, 0 ≤ a → a ≤ b → g a ≤ g b := by
    intro a b ha hab
    have hba : (0:ℝ) ≤ b - a := by linarith
    have h := hadd a (b - a) ha hba
    have hb : a + (b - a) = b := by ring
    rw [hb] at h
    have := hnn (b - a) hba
    linarith
  -- natuurlijke schaling
  have hnat : ∀ n : ℕ, g n = n * g 1 := by
    intro n
    simpa using iterate_add g hadd 1 (by norm_num) n
  -- rationale schaling
  have hrat : ∀ q : ℚ, 0 ≤ q → g q = q * g 1 := by
    intro q hq
    have hqR : (0:ℝ) ≤ (q:ℝ) := by exact_mod_cast hq
    have hden_pos : (0:ℝ) < (q.den : ℝ) := by exact_mod_cast q.den_pos
    have hden_ne : ((q.den : ℝ)) ≠ 0 := ne_of_gt hden_pos
    have hnum0 : 0 ≤ q.num := Rat.num_nonneg.mpr hq
    have h1 : g ((q.den : ℝ) * (q:ℝ)) = (q.den : ℕ) * g (q:ℝ) := by
      simpa using iterate_add g hadd (q:ℝ) hqR q.den
    have h2 : ((q.den : ℝ)) * (q:ℝ) = (q.num : ℝ) := by
      rw [mul_comm]
      exact_mod_cast Rat.mul_den_eq_num q
    have h3 : g ((q.num : ℝ)) = (q.num : ℝ) * g 1 := by
      have hn : ((q.num.toNat : ℕ) : ℝ) = (q.num : ℝ) := by
        exact_mod_cast Int.toNat_of_nonneg hnum0
      rw [← hn]
      simpa using hnat q.num.toNat
    have hkey : (q.den : ℝ) * g (q:ℝ) = (q.den : ℝ) * ((q:ℝ) * g 1) := by
      have hL : (q.den : ℝ) * g (q:ℝ) = (q.num : ℝ) * g 1 := by
        rw [← h1, h2, h3]
      have hR : (q.den : ℝ) * ((q:ℝ) * g 1) = (q.num : ℝ) * g 1 := by
        rw [← mul_assoc, h2]
      rw [hL, hR]
    exact mul_left_cancel₀ hden_ne hkey
  -- sandwich
  intro s hs
  rcases eq_or_lt_of_le (hnn 1 (by norm_num)) with h10 | h1pos
  · -- g 1 = 0 (h10 : 0 = g 1) ⇒ g s = 0
    have hub : g s ≤ g ((⌈s⌉₊ : ℕ) : ℝ) := hmono s _ hs (Nat.le_ceil s)
    have hceil : g ((⌈s⌉₊ : ℕ) : ℝ) = 0 := by rw [hnat ⌈s⌉₊, ← h10]; ring
    have hlb := hnn s hs
    have hgs : g s = 0 := le_antisymm (by rw [← hceil]; exact hub) hlb
    rw [hgs, ← h10]; ring
  · apply le_antisymm
    · -- g s ≤ g 1 * s (+ε)
      apply le_of_forall_pos_le_add
      intro ε hε
      have hlt : s < s + ε / g 1 := by
        have := div_pos hε h1pos; linarith
      obtain ⟨r, hr1, hr2⟩ := exists_rat_btwn hlt
      have hr0 : (0:ℚ) ≤ r := by
        have : (0:ℝ) ≤ (r:ℝ) := le_trans hs (le_of_lt hr1)
        exact_mod_cast this
      have hmr : g s ≤ (r:ℝ) * g 1 := by
        have h := hmono s r hs (le_of_lt hr1)
        rwa [hrat r hr0] at h
      have hgr : (r:ℝ) * g 1 < (s + ε / g 1) * g 1 :=
        mul_lt_mul_of_pos_right hr2 h1pos
      have hexp : (s + ε / g 1) * g 1 = g 1 * s + ε := by
        field_simp
      linarith
    · -- g 1 * s ≤ g s (+ε)
      apply le_of_forall_pos_le_add
      intro ε hε
      by_cases hcase : ε / g 1 ≤ s
      · have hlt : s - ε / g 1 < s := by
          have := div_pos hε h1pos; linarith
        obtain ⟨q, hq1, hq2⟩ := exists_rat_btwn hlt
        have hq0 : (0:ℚ) ≤ q := by
          have h0 : (0:ℝ) ≤ s - ε / g 1 := by linarith
          have : (0:ℝ) ≤ (q:ℝ) := le_trans h0 (le_of_lt hq1)
          exact_mod_cast this
        have hqs : (q:ℝ) * g 1 ≤ g s := by
          have h := hmono q s (by exact_mod_cast hq0) (le_of_lt hq2)
          rwa [hrat q hq0] at h
        have hlow : (s - ε / g 1) * g 1 < (q:ℝ) * g 1 :=
          mul_lt_mul_of_pos_right hq1 h1pos
        have hexp : (s - ε / g 1) * g 1 = g 1 * s - ε := by
          field_simp
        linarith
      · have hcase' : s < ε / g 1 := not_le.mp hcase
        have h1 : g 1 * s < ε := by
          have h2 := mul_lt_mul_of_pos_left hcase' h1pos
          have h3 : g 1 * (ε / g 1) = ε := by field_simp
          linarith
        linarith [hnn s hs]

/-- De EEML-kernstelling: compositie + codomein ⇒ P = exp(−λs). -/
theorem eeml_kern (P : ℝ → ℝ)
    (hpos : ∀ s : ℝ, 0 ≤ s → 0 < P s)
    (hle  : ∀ s : ℝ, 0 ≤ s → P s ≤ 1)
    (hmul : ∀ s t : ℝ, 0 ≤ s → 0 ≤ t → P (s + t) = P s * P t) :
    ∃ lam : ℝ, 0 ≤ lam ∧ ∀ s : ℝ, 0 ≤ s → P s = Real.exp (-(lam * s)) := by
  set g : ℝ → ℝ := fun s => -Real.log (P s) with hg
  have hadd : ∀ s t : ℝ, 0 ≤ s → 0 ≤ t → g (s + t) = g s + g t := by
    intro s t hs ht
    simp only [hg]
    rw [hmul s t hs ht, Real.log_mul (ne_of_gt (hpos s hs)) (ne_of_gt (hpos t ht))]
    ring
  have hnn : ∀ s : ℝ, 0 ≤ s → 0 ≤ g s := by
    intro s hs
    simp only [hg, neg_nonneg]
    exact Real.log_nonpos (le_of_lt (hpos s hs)) (hle s hs)
  have hlin := additive_nonneg_linear g hadd hnn
  refine ⟨g 1, hnn 1 (by norm_num), ?_⟩
  intro s hs
  have hgl := hlin s hs
  have hlog : Real.log (P s) = -(g 1 * s) := by
    have : g s = -Real.log (P s) := by simp only [hg]
    linarith [hgl, this]
  calc P s = Real.exp (Real.log (P s)) := (Real.exp_log (hpos s hs)).symm
    _ = Real.exp (-(g 1 * s)) := by rw [hlog]

#print axioms EEML.eeml_kern

end EEML