MaxLeafFormal.CCRate

Source file MaxLeafFormal/CCRate.lean of the Lean project (how it is built).

Statements of the paper formalised here: Lemma 11.12, Lemma D.9.8.

import MaxLeafFormal.CCData
import MaxLeafFormal.EncloseH

/-!
# Graph Motif: the rate certificate for `c = 2`, and the checks for `c = 3`

`lem:cc-cert` of the paper: every
`σ ∈ (0, 1]` satisfies (R1), (R2) or (R3) of `cor:cc-rate` with `c = 2`, `ω = 0.71755`, `Γ₀ = 16`,
and `0.71755 < log₂ 1.6444` (`Real.logb 2`). For `c = 3`, `cc_c3` proves the arithmetic facts that
`ancillary/component_count_rate_interval.py` checks: `spr(6/5) = 3`, `(6/5)/3 < 3`, `1 + 2·(6/5)/3 = 9/5`,
`3^64 < 5^45` and `log₂(9/5) < H(26/45)`. It does not state `cor:cc-T`(b) itself; the manuscript
proves (B2) and (B3) by hand.
-/

namespace MaxLeafFormal.CCRate

open Enclose

/-- `spr(M_ζ) = (1 + 2ζ + √(1 + 4ζ²))/2`. -/
noncomputable def spr (ζ : ℝ) : ℝ := (1 + 2 * ζ + Real.sqrt (1 + 4 * ζ^2)) / 2

/-- `φ_{c,σ}(y) = log(1+y) - σ log(spr(M_{cy/2})/c)` (`eq:cc-phi`). -/
noncomputable def phi (c σ y : ℝ) : ℝ := Real.logb 2 (1 + y) - σ * Real.logb 2 (spr (c * y / 2) / c)

def R1 (c ω σ : ℝ) : Prop := σ * Real.logb 2 c ≤ ω

def R2 (ω σ : ℝ) : Prop := MLS.H2 σ ≤ ω

def R3 (c ω Γ₀ σ : ℝ) : Prop :=
  ∃ y > 0, phi c σ y ≤ ω ∧ spr (c * y / 2) / min 1 (c * y / 2) ≤ Γ₀

/-! ### The checker (over `ℚ`) -/

/-- An enclosure of `spr(ζ)` for rational `ζ ≥ 0`. -/
def sprQ (ζ : ℚ) : QI := ((QI.pt (1 + 2 * ζ)).add (sqrtQ (1 + 4 * ζ ^ 2) 80)).smul (1 / 2)

/-- The interval `[a, a + 10⁻⁴]` with parameter `y` (`c = 2`, so `ζ = y`): `y > 0`,
`max(φ_{2,a}(y), φ_{2,a+10⁻⁴}(y)) < 0.71755` with enclosures, and `spr(y)/min(1, y) ≤ 16`. -/
def intervalOK (a y : ℚ) : Bool :=
  let b := a + 1 / 10 ^ 4
  let S := sprQ y
  let L1 := log2Q (1 + y) 60 80
  let LS := log2I (S.smul (1 / 2)) 60 80
  decide (0 < y) && decide (0 < S.lo) &&
    decide (L1.hi - a * LS.lo < 71755 / 10 ^ 5) && decide (L1.hi - b * LS.lo < 71755 / 10 ^ 5) &&
    decide (0 ≤ a) && decide (S.hi ≤ 16 * min 1 y)

/-- Consecutive intervals of width `10⁻⁴` starting at `a`. -/
def checkIntervals : ℚ → List ℚ → Bool
  | _, [] => true
  | a, y :: rest => intervalOK a y && checkIntervals (a + 1 / 10 ^ 4) rest

/-! The kernel check, in chunks of 100 intervals (chunk `k` starts at `0.7175 + 100k/10^4`). -/
set_option maxRecDepth 100000 in

theorem intervals_ok_0 :
    checkIntervals (7175 / 10 ^ 4 + 0 / 10 ^ 4) ((ys.drop 0).take 100) = true := by
  decide +kernel
set_option maxRecDepth 100000 in

theorem intervals_ok_1 :
    checkIntervals (7175 / 10 ^ 4 + 100 / 10 ^ 4) ((ys.drop 100).take 100) = true := by
  decide +kernel
set_option maxRecDepth 100000 in

theorem intervals_ok_2 :
    checkIntervals (7175 / 10 ^ 4 + 200 / 10 ^ 4) ((ys.drop 200).take 100) = true := by
  decide +kernel
set_option maxRecDepth 100000 in

theorem intervals_ok_3 :
    checkIntervals (7175 / 10 ^ 4 + 300 / 10 ^ 4) ((ys.drop 300).take 100) = true := by
  decide +kernel
set_option maxRecDepth 100000 in

theorem intervals_ok_4 :
    checkIntervals (7175 / 10 ^ 4 + 400 / 10 ^ 4) ((ys.drop 400).take 100) = true := by
  decide +kernel
set_option maxRecDepth 100000 in

theorem intervals_ok_5 :
    checkIntervals (7175 / 10 ^ 4 + 500 / 10 ^ 4) ((ys.drop 500).take 100) = true := by
  decide +kernel
set_option maxRecDepth 100000 in

theorem intervals_ok_6 :
    checkIntervals (7175 / 10 ^ 4 + 600 / 10 ^ 4) ((ys.drop 600).take 100) = true := by
  decide +kernel
set_option maxRecDepth 100000 in

theorem intervals_ok_7 :
    checkIntervals (7175 / 10 ^ 4 + 700 / 10 ^ 4) ((ys.drop 700).take 100) = true := by
  decide +kernel
set_option maxRecDepth 100000 in

theorem intervals_ok_8 :
    checkIntervals (7175 / 10 ^ 4 + 800 / 10 ^ 4) ((ys.drop 800).take 100) = true := by
  decide +kernel
set_option maxRecDepth 100000 in

theorem intervals_ok_9 :
    checkIntervals (7175 / 10 ^ 4 + 900 / 10 ^ 4) ((ys.drop 900).take 100) = true := by
  decide +kernel

theorem ys_length : ys.length = 925 := by decide +kernel

/-! ### Soundness and the theorems -/

/-- Every entry of a checked list is an accepted interval. -/
theorem checkIntervals_sound (l : List ℚ) : ∀ (a : ℚ), checkIntervals a l = true →
    ∀ i, i < l.length → intervalOK (a + (i : ℚ) / 10 ^ 4) (l.getD i 0) = true := by
  induction l with
  | nil => intro a _ i hi; exact absurd hi (Nat.not_lt_zero _)
  | cons y rest ih =>
    intro a h i hi
    simp only [checkIntervals, Bool.and_eq_true] at h
    cases i with
    | zero => simpa using h.1
    | succ i =>
      have ih' := ih (a + 1 / 10 ^ 4) h.2 i (by simpa using hi)
      have e : a + ((i + 1 : ℕ) : ℚ) / 10 ^ 4 = a + 1 / 10 ^ 4 + (i : ℚ) / 10 ^ 4 := by
        push_cast; ring
      rw [e]
      simpa using ih'

/-- A chunk `s` covers the indices `s ≤ j < s + 100`. -/
theorem chunk (s : ℕ) (a : ℚ) (ha : a = 7175 / 10 ^ 4 + (s : ℚ) / 10 ^ 4)
    (hs : checkIntervals a ((ys.drop s).take 100) = true) (j : ℕ)
    (h1 : s ≤ j) (h2 : j < s + 100) (hj : j < 925) :
    intervalOK (7175 / 10 ^ 4 + (j : ℚ) / 10 ^ 4) (ys.getD j 0) = true := by
  have hlen : j - s < ((ys.drop s).take 100).length := by
    simp only [List.length_take, List.length_drop, ys_length]; omega
  have key := checkIntervals_sound _ a hs (j - s) hlen
  have e1 : a + ((j - s : ℕ) : ℚ) / 10 ^ 4 = 7175 / 10 ^ 4 + (j : ℚ) / 10 ^ 4 := by
    rw [ha, Nat.cast_sub h1]; ring
  have e2 : ((ys.drop s).take 100).getD (j - s) 0 = ys.getD j 0 := by
    rw [List.getD_eq_getElem?_getD, List.getD_eq_getElem?_getD, List.getElem?_take,
      ite_eq_left (show j - s < 100 by omega), List.getElem?_drop, show s + (j - s) = j by omega]
  rw [e1, e2] at key
  exact key

theorem all_ok (j : ℕ) (hj : j < 925) :
    intervalOK (7175 / 10 ^ 4 + (j : ℚ) / 10 ^ 4) (ys.getD j 0) = true := by
  obtain ⟨k, hk⟩ : ∃ k, j / 100 = k := ⟨_, rfl⟩
  have hk10 : k < 10 := by omega
  interval_cases k
  · exact chunk 0 _ (by norm_num) intervals_ok_0 j (by omega) (by omega) hj
  · exact chunk 100 _ (by norm_num) intervals_ok_1 j (by omega) (by omega) hj
  · exact chunk 200 _ (by norm_num) intervals_ok_2 j (by omega) (by omega) hj
  · exact chunk 300 _ (by norm_num) intervals_ok_3 j (by omega) (by omega) hj
  · exact chunk 400 _ (by norm_num) intervals_ok_4 j (by omega) (by omega) hj
  · exact chunk 500 _ (by norm_num) intervals_ok_5 j (by omega) (by omega) hj
  · exact chunk 600 _ (by norm_num) intervals_ok_6 j (by omega) (by omega) hj
  · exact chunk 700 _ (by norm_num) intervals_ok_7 j (by omega) (by omega) hj
  · exact chunk 800 _ (by norm_num) intervals_ok_8 j (by omega) (by omega) hj
  · exact chunk 900 _ (by norm_num) intervals_ok_9 j (by omega) (by omega) hj

theorem sprQ_sound (ζ : ℚ) : (sprQ ζ).Contains (spr (ζ : ℝ)) := by
  have h := QI.smul_sound (1 / 2 : ℚ)
    (QI.add_sound (QI.pt_sound (1 + 2 * ζ)) (sqrtQ_sound (q := 1 + 4 * ζ ^ 2) (by positivity) 80))
  have e : spr (ζ : ℝ) = ((1 / 2 : ℚ) : ℝ) *
      (((1 + 2 * ζ : ℚ) : ℝ) + Real.sqrt ((1 + 4 * ζ ^ 2 : ℚ) : ℝ)) := by
    unfold spr; push_cast; ring
  rw [e]
  exact h

/-- The real core of `R3` for one interval: `σ ∈ [a, b]`, affine in `σ`. -/
theorem R3_real {σ a b y L1 LS : ℝ} (hσ0 : 0 ≤ σ) (hσ1 : a ≤ σ) (hσ2 : σ ≤ b)
    (hA : Real.logb 2 (1 + y) ≤ L1) (hB : LS ≤ Real.logb 2 (spr y / 2))
    (hC : L1 - a * LS < 71755 / 10 ^ 5) (hD : L1 - b * LS < 71755 / 10 ^ 5) :
    Real.logb 2 (1 + y) - σ * Real.logb 2 (spr y / 2) ≤ 71755 / 10 ^ 5 := by
  have h2 := mul_le_mul_of_nonneg_left hB hσ0
  rcases le_total 0 LS with hL | hL
  · have h1 := mul_le_mul_of_nonneg_right hσ1 hL
    linarith
  · have h1 := mul_le_mul_of_nonpos_right hσ2 hL
    linarith

theorem intervalOK_R3 {a y : ℚ} (h : intervalOK a y = true) {σ : ℝ}
    (hσ1 : (a : ℝ) ≤ σ) (hσ2 : σ ≤ ((a + 1 / 10 ^ 4 : ℚ) : ℝ)) :
    R3 2 (71755 / 10 ^ 5) 16 σ := by
  simp only [intervalOK, Bool.and_eq_true, decide_eq_true_eq] at h
  obtain ⟨⟨⟨⟨⟨hy, hS0⟩, hL2⟩, hL3⟩, ha0⟩, hS1⟩ := h
  have hc : ((sprQ y).smul (1 / 2)).lo = 1 / 2 * (sprQ y).lo := by
    unfold QI.smul; rw [ite_eq_left (show (0:ℚ) ≤ 1 / 2 by norm_num)]
  have hpos : 0 < ((sprQ y).smul (1 / 2)).lo := by
    rw [hc]; exact mul_pos (by norm_num) hS0
  have hLS := log2I_sound (QI.smul_sound (1 / 2 : ℚ) (sprQ_sound y)) hpos 60 80
  have hA := log2Q_sound (q := 1 + y) (by linarith) 60 80
  have hS := sprQ_sound y
  simp only [QI.Contains] at hLS hA hS
  have e2 : ((1 / 2 : ℚ) : ℝ) * spr (y : ℝ) = spr (y : ℝ) / 2 := by push_cast; ring
  rw [e2] at hLS
  have hA' : Real.logb 2 (1 + (y : ℝ)) ≤ ((log2Q (1 + y) 60 80).hi : ℝ) := by
    have h := hA.2; push_cast at h; exact h
  have hσ2' : σ ≤ (a : ℝ) + 1 / 10 ^ 4 := by push_cast at hσ2; exact hσ2
  have hC : ((log2Q (1 + y) 60 80).hi : ℝ) - (a : ℝ) *
      ((log2I ((sprQ y).smul (1 / 2)) 60 80).lo : ℝ) < 71755 / 10 ^ 5 := by
    have := (Rat.cast_lt (K := ℝ)).mpr hL2
    push_cast at this; linarith
  have hD : ((log2Q (1 + y) 60 80).hi : ℝ) - ((a : ℝ) + 1 / 10 ^ 4) *
      ((log2I ((sprQ y).smul (1 / 2)) 60 80).lo : ℝ) < 71755 / 10 ^ 5 := by
    have := (Rat.cast_lt (K := ℝ)).mpr hL3
    push_cast at this; linarith
  have hσ0 : (0 : ℝ) ≤ σ := le_trans (Rat.cast_nonneg.mpr ha0) hσ1
  unfold R3
  refine ⟨(y : ℝ), by exact_mod_cast hy, ?_, ?_⟩
  · unfold phi
    rw [show (2 : ℝ) * (y : ℝ) / 2 = (y : ℝ) by ring]
    exact R3_real hσ0 hσ1 hσ2' hA' hLS.1 hC hD
  · rw [show (2 : ℝ) * (y : ℝ) / 2 = (y : ℝ) by ring]
    have hm : (0 : ℝ) < min 1 (y : ℝ) := lt_min one_pos (by exact_mod_cast hy)
    have hS1' : ((sprQ y).hi : ℝ) ≤ 16 * min 1 (y : ℝ) := by
      have := (Rat.cast_le (K := ℝ)).mpr hS1
      push_cast at this; exact this
    rw [div_le_iff₀ hm]
    linarith [hS.2]

/-- `spr(M_ζ)` is a root of `λ² - (1 + 2ζ)λ + ζ` (the characteristic polynomial of `M_ζ`). -/
theorem spr_root (ζ : ℝ) : spr ζ ^ 2 - (1 + 2 * ζ) * spr ζ + ζ = 0 := by
  have hs : Real.sqrt (1 + 4 * ζ ^ 2) ^ 2 = 1 + 4 * ζ ^ 2 :=
    Real.sq_sqrt (by positivity)
  unfold spr
  linear_combination (1 / 4 : ℝ) * hs

theorem cc_cert :
    (∀ σ ∈ Set.Ioc (0:ℝ) 1, R1 2 (71755/10^5) σ ∨ R2 (71755/10^5) σ ∨ R3 2 (71755/10^5) 16 σ) ∧
    (71755/10^5 : ℝ) < Real.logb 2 (16444/10^4) := by
  refine ⟨?_, ?_⟩
  · intro σ hσ
    obtain ⟨hσ0, hσ1⟩ := Set.mem_Ioc.mp hσ
    rcases le_total σ (7175 / 10 ^ 4) with h1 | h1
    · left
      unfold R1
      rw [Real.logb_self_eq_one (by norm_num)]
      linarith
    rcases lt_or_ge σ (81 / 100) with h2 | h2
    · right; right
      have ht0 : 0 ≤ (σ - 7175 / 10 ^ 4) * 10 ^ 4 := mul_nonneg (by linarith) (by norm_num)
      obtain ⟨j, hj⟩ : ∃ j : ℕ, j = ⌊(σ - 7175 / 10 ^ 4) * 10 ^ 4⌋₊ := ⟨_, rfl⟩
      have hf1 : (j : ℝ) ≤ (σ - 7175 / 10 ^ 4) * 10 ^ 4 := by
        rw [hj]; exact Nat.floor_le ht0
      have hf2 : (σ - 7175 / 10 ^ 4) * 10 ^ 4 < (j : ℝ) + 1 := by
        rw [hj]; exact Nat.lt_floor_add_one _
      have hj925 : j < 925 := by
        have : (j : ℝ) < 925 := by linarith
        exact_mod_cast this
      have hσa : ((7175 / 10 ^ 4 + (j : ℚ) / 10 ^ 4 : ℚ) : ℝ) ≤ σ := by
        push_cast; linarith
      have hσb : σ ≤ ((7175 / 10 ^ 4 + (j : ℚ) / 10 ^ 4 + 1 / 10 ^ 4 : ℚ) : ℝ) := by
        push_cast; linarith
      exact intervalOK_R3 (all_ok j hj925) hσa hσb
    · right; left
      unfold R2
      have hanti : Real.binEntropy σ ≤ Real.binEntropy (81 / 100) :=
        Real.binEntropy_strictAntiOn.antitoneOn
          ⟨by norm_num, by norm_num⟩ ⟨le_trans (by norm_num) h2, hσ1⟩ h2
      have hH : MLS.H2 σ ≤ MLS.H2 (81 / 100) := by
        unfold MLS.H2
        exact div_le_div_of_nonneg_right hanti (Real.log_nonneg (by norm_num))
      have hQ := H2Q_sound (x := 81 / 100) (by norm_num) (by norm_num) 60 80
      have hb : (H2Q (81 / 100) 60 80).hi < 71755 / 10 ^ 5 := by decide +kernel
      have hc : ((81 / 100 : ℚ) : ℝ) = 81 / 100 := by norm_num
      simp only [QI.Contains] at hQ
      rw [hc] at hQ
      have hbR := (Rat.cast_lt (K := ℝ)).mpr hb
      push_cast at hbR
      linarith [hQ.2]
  · have hL := log2Q_sound (q := 16444 / 10 ^ 4) (by norm_num) 60 80
    have hlo : (71755 / 10 ^ 5 : ℚ) < (log2Q (16444 / 10 ^ 4) 60 80).lo := by
      decide +kernel
    have hlo' := (Rat.cast_lt (K := ℝ)).mpr hlo
    push_cast at hlo'
    have hc : ((16444 / 10 ^ 4 : ℚ) : ℝ) = 16444 / 10 ^ 4 := by norm_num
    simp only [QI.Contains] at hL
    rw [hc] at hL
    linarith [hL.1]

theorem cc_c3 :
    spr (6/5) = 3 ∧ (6/5 : ℝ) / 3 < 3 ∧ (1 : ℝ) + 2 * (6/5) / 3 = 9/5 ∧
    (3:ℕ) ^ 64 < 5 ^ 45 ∧ Real.logb 2 (9/5) < MLS.H2 (26/45) := by
  have hsq : Real.sqrt (1 + 4 * (6 / 5 : ℝ) ^ 2) = 13 / 5 := by
    rw [show (1 : ℝ) + 4 * (6 / 5) ^ 2 = (13 / 5) ^ 2 by norm_num]
    exact Real.sqrt_sq (by norm_num)
  have hspr : spr (6 / 5) = 3 := by
    unfold spr; rw [hsq]; norm_num
  refine ⟨hspr, by norm_num, by norm_num, by norm_num, ?_⟩
  have hL := log2Q_sound (q := 9 / 5) (by norm_num) 60 80
  have hH := H2Q_sound (x := 26 / 45) (by norm_num) (by norm_num) 60 80
  have hcmp : (log2Q (9 / 5) 60 80).hi < (H2Q (26 / 45) 60 80).lo := by decide +kernel
  have hcmpR : (((log2Q (9 / 5) 60 80).hi : ℚ) : ℝ) < (((H2Q (26 / 45) 60 80).lo : ℚ) : ℝ) := by
    exact_mod_cast hcmp
  have e1 : ((9 / 5 : ℚ) : ℝ) = 9 / 5 := by norm_num
  have e2 : ((26 / 45 : ℚ) : ℝ) = 26 / 45 := by norm_num
  simp only [QI.Contains] at hL hH
  rw [e1] at hL
  rw [e2] at hH
  linarith [hL.2, hH.1]

end MaxLeafFormal.CCRate