MaxLeafFormal.PDRate

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

Statements of the paper formalised here: Lemma D.11.19, Theorem 12.10.

import MaxLeafFormal.Enclose

/-!
# Paired Domination: the rate

`lem:fa-pd-rate`, first claim, of the paper: `max_{n̄ ∈ {2,3,…}} h_{n̄}(y*) = h_6(y*) = ω_pr = 0.7050097…`, and the constants of
`thm:fa-pd`: `2^{ω_pr} = 1.6301556…`, `ω_pr < 0.70502`, `2^{0.70502} < 1.6302`.
-/

namespace MaxLeafFormal.PDRate

open Enclose

/-- `h_{n̄}(y) = log(1+y) - (1/n̄) log(y² + 4·2^{-n̄/2} y + 2^{-n̄})`. -/
noncomputable def hpd (n y : ℝ) : ℝ :=
  Real.logb 2 (1 + y) - (1 / n) * Real.logb 2 (y^2 + 4 * (2:ℝ) ^ (-n / 2) * y + (2:ℝ) ^ (-n))

/-- `y* = (3√3 - 1)/16`. -/
noncomputable def ystar : ℝ := (3 * Real.sqrt 3 - 1) / 16

/-- `ω_pr = h_6(y*)`. -/
noncomputable def ωpr : ℝ := hpd 6 ystar

/-! ### Rational enclosures (helpers) -/

/-- Enclosure of `y* = (3√3 - 1)/16`. -/
def pdY : QI := QI.smul (1 / 16) (QI.add (QI.smul 3 (sqrtQ 3 80)) (QI.pt (-1)))

/-- Enclosure of `y² + 4·2^{-m/2}·y + 2^{-m}` for `y ∈ Y`. -/
def pdArg (m : ℕ) (Y : QI) : QI :=
  QI.add (QI.add (QI.mul Y Y) (QI.mul (QI.smul 4 (pow2Q (-(m:ℚ) / 2) 60 80)) Y))
    (pow2Q (-(m:ℚ)) 60 80)

/-- Enclosure of `h_m(y)` for `y ∈ Y`. -/
def pdH (m : ℕ) (Y : QI) : QI :=
  QI.sub (log2I (QI.add (QI.pt 1) Y) 60 80)
    (QI.smul (1 / (m:ℚ)) (log2I (pdArg m Y) 60 80))

/-- Enclosure of `log₂(1+y*)`-part and `log₂ y*`-part of the tail bound. -/
def pdTail : ℚ :=
  (log2I (QI.add (QI.pt 1) pdY) 60 80).hi + 2 / 41 * (-(log2I pdY 60 80).lo)

theorem pdY_sound : pdY.Contains ystar := by
  have hS := sqrtQ_sound (q := 3) (by norm_num) 80
  have h3 := QI.smul_sound 3 hS
  have h4 := QI.add_sound h3 (QI.pt_sound (-1))
  have h5 := QI.smul_sound (1 / 16) h4
  have e : ystar = ((1 / 16 : ℚ) : ℝ) * (((3 : ℚ) : ℝ) * Real.sqrt ((3 : ℚ) : ℝ) + ((-1 : ℚ) : ℝ)) := by
    unfold ystar; push_cast <;> ring
  rw [e]
  exact h5

/-- Rational order passes to the reals (lower bound with a literal). -/
theorem qR_le_of {a c : ℚ} {b : ℝ} (h : a ≤ c) (hb : (a:ℝ) = b) : b ≤ (c:ℝ) :=
  hb ▸ Rat.cast_le.mpr h

/-- Rational order passes to the reals (upper bound with a literal). -/
theorem qR_lt_of {a c : ℚ} {d : ℝ} (h : a < c) (hd : (c:ℝ) = d) : (a:ℝ) < d :=
  hd ▸ Rat.cast_lt.mpr h

theorem pdY_lo_pos : (0:ℚ) < pdY.lo := by decide +kernel

theorem pdY_hi_lt : pdY.hi < 1 := by decide +kernel

theorem pdArg_sound {Y : QI} {y : ℝ} (hY : Y.Contains y) (m : ℕ) :
    (pdArg m Y).Contains (y ^ 2 + 4 * (2:ℝ) ^ (-(m:ℝ) / 2) * y + (2:ℝ) ^ (-(m:ℝ))) := by
  have hA1 := QI.mul_sound hY hY
  have hP := pow2Q_sound (-(m:ℚ) / 2) 60 80
  have hc1 : ((-(m:ℚ) / 2 : ℚ) : ℝ) = -(m:ℝ) / 2 := by push_cast <;> ring
  rw [hc1] at hP
  have hS := QI.smul_sound 4 hP
  have hM := QI.mul_sound hS hY
  have hQ := pow2Q_sound (-(m:ℚ)) 60 80
  have hc2 : ((-(m:ℚ) : ℚ) : ℝ) = -(m:ℝ) := by push_cast <;> ring
  rw [hc2] at hQ
  have hsum := QI.add_sound (QI.add_sound hA1 hM) hQ
  have e : y ^ 2 + 4 * (2:ℝ) ^ (-(m:ℝ) / 2) * y + (2:ℝ) ^ (-(m:ℝ))
      = y * y + ((4 : ℚ) : ℝ) * (2:ℝ) ^ (-(m:ℝ) / 2) * y + (2:ℝ) ^ (-(m:ℝ)) := by
    push_cast <;> ring
  rw [e]
  exact hsum

theorem pdH_sound {m : ℕ} {Y : QI} {y : ℝ} (hY : Y.Contains y) (hY0 : 0 < Y.lo)
    (hA : 0 < (pdArg m Y).lo) : (pdH m Y).Contains (hpd m y) := by
  have hJ : (QI.add (QI.pt 1) Y).Contains (1 + y) := by
    have := QI.add_sound (QI.pt_sound 1) hY
    push_cast at this
    exact this
  have hJ0 : 0 < (QI.add (QI.pt 1) Y).lo := by
    show (0:ℚ) < 1 + Y.lo
    linarith
  have hL1 := log2I_sound hJ hJ0 60 80
  have hL2 := log2I_sound (pdArg_sound hY m) hA 60 80
  have hS := QI.smul_sound (1 / (m:ℚ)) hL2
  have hD := QI.sub_sound hL1 hS
  have hc : hpd m y = Real.logb 2 (1 + y) - ((1 / (m:ℚ) : ℚ) : ℝ) *
      Real.logb 2 (y ^ 2 + 4 * (2:ℝ) ^ (-(m:ℝ) / 2) * y + (2:ℝ) ^ (-(m:ℝ))) := by
    unfold hpd; push_cast <;> ring
  rw [hc]
  exact hD

/-- Tail bound: for `m ≥ 41` and `0 < y ≤ 1`. -/
theorem hpd_tail {m : ℕ} (hm : 41 ≤ m) {y : ℝ} (hy0 : 0 < y) (hy1 : y ≤ 1) :
    hpd m y ≤ Real.logb 2 (1 + y) + 2 / 41 * (-(Real.logb 2 y)) := by
  have hmr : (41:ℝ) ≤ (m:ℝ) := by exact_mod_cast hm
  have h1 : (0:ℝ) ≤ (2:ℝ) ^ (-(m:ℝ) / 2) := by positivity
  have h2 : (0:ℝ) ≤ (2:ℝ) ^ (-(m:ℝ)) := by positivity
  have h3 : 0 ≤ 4 * (2:ℝ) ^ (-(m:ℝ) / 2) * y :=
    mul_nonneg (mul_nonneg (by norm_num) h1) hy0.le
  have hA : y ^ 2 ≤ y ^ 2 + 4 * (2:ℝ) ^ (-(m:ℝ) / 2) * y + (2:ℝ) ^ (-(m:ℝ)) := by
    linarith
  have hpos : 0 < y ^ 2 := pow_pos hy0 2
  have hlogsq : Real.logb 2 (y ^ 2) = 2 * Real.logb 2 y := by
    unfold Real.logb
    rw [Real.log_pow]
    push_cast <;> ring
  have hLA : 2 * Real.logb 2 y ≤
      Real.logb 2 (y ^ 2 + 4 * (2:ℝ) ^ (-(m:ℝ) / 2) * y + (2:ℝ) ^ (-(m:ℝ))) := by
    rw [← hlogsq]
    exact Real.logb_le_logb_of_le (by norm_num) hpos hA
  have hL : Real.logb 2 y ≤ 0 := Real.logb_nonpos (by norm_num) hy0.le hy1
  have hinv : 1 / (m:ℝ) ≤ 1 / 41 := one_div_le_one_div_of_le (by norm_num) hmr
  have hm1 : 0 ≤ 1 / (m:ℝ) := by positivity
  have hnn : 0 ≤ -(Real.logb 2 y) := by linarith
  have e1 := mul_le_mul_of_nonneg_left hLA hm1
  have e2 := mul_le_mul_of_nonneg_right hinv hnn
  unfold hpd
  linarith [e1, e2]

theorem finite_check : ∀ m ∈ List.range 41, 2 ≤ m → m ≠ 6 →
    0 < (pdArg m pdY).lo ∧ (pdH m pdY).hi < (pdH 6 pdY).lo := by
  decide +kernel

theorem pdA6_pos : 0 < (pdArg 6 pdY).lo := by decide +kernel

theorem pdH6_lo : (7050097/10^7 : ℚ) ≤ (pdH 6 pdY).lo := by decide +kernel

theorem pdH6_hi : (pdH 6 pdY).hi < 7050098/10^7 := by decide +kernel

theorem pdTail_lt : pdTail < 7050097/10^7 := by decide +kernel

theorem pdH6_contains : (pdH 6 pdY).Contains ωpr := by
  have h := pdH_sound (m := 6) pdY_sound pdY_lo_pos pdA6_pos
  rw [show ((6:ℕ):ℝ) = (6:ℝ) by norm_num] at h
  exact h

theorem pdC1 : (16301556/10^7 : ℚ) ≤ (pow2Q (pdH 6 pdY).lo 60 80).lo := by decide +kernel

theorem pdC2 : (pow2Q (pdH 6 pdY).hi 60 80).hi < 16301557/10^7 := by decide +kernel

theorem pdC3 : (pdH 6 pdY).hi < 70502/10^5 := by decide +kernel

theorem pdC4 : (pow2Q (70502/10^5) 60 80).hi < 16302/10^4 := by decide +kernel

theorem pd_rate : (∀ m : ℕ, 2 ≤ m → hpd m ystar ≤ ωpr) ∧
    (7050097/10^7 : ℝ) ≤ ωpr ∧ ωpr < 7050098/10^7 := by
  have hW' : ((pdH 6 pdY).lo : ℝ) ≤ ωpr ∧ ωpr ≤ ((pdH 6 pdY).hi : ℝ) := pdH6_contains
  have c1 : (7050097/10^7 : ℝ) ≤ ((pdH 6 pdY).lo : ℝ) := qR_le_of pdH6_lo (by norm_num)
  have c2 : ((pdH 6 pdY).hi : ℝ) < 7050098/10^7 := qR_lt_of pdH6_hi (by norm_num)
  have hWl : (7050097/10^7 : ℝ) ≤ ωpr := le_trans c1 hW'.1
  have hWu : ωpr < 7050098/10^7 := lt_of_le_of_lt hW'.2 c2
  refine ⟨?_, hWl, hWu⟩
  intro m hm
  by_cases h41 : 41 ≤ m
  · have hYs := pdY_sound
    have hY0r : (0:ℝ) < ((pdY.lo : ℚ) : ℝ) := by exact_mod_cast pdY_lo_pos
    have hY1r : ((pdY.hi : ℚ) : ℝ) < 1 := by exact_mod_cast pdY_hi_lt
    have hypos : 0 < ystar := lt_of_lt_of_le hY0r hYs.1
    have hyle : ystar ≤ 1 := le_trans hYs.2 hY1r.le
    have ht := hpd_tail h41 hypos hyle
    have hI : (QI.add (QI.pt 1) pdY).Contains (1 + ystar) := by
      have := QI.add_sound (QI.pt_sound 1) hYs
      push_cast at this
      exact this
    have hJ0 : 0 < (QI.add (QI.pt 1) pdY).lo := by
      show (0:ℚ) < 1 + pdY.lo
      linarith [pdY_lo_pos]
    have hb1 := (log2I_sound hI hJ0 60 80).2
    have hb2 := (log2I_sound hYs pdY_lo_pos 60 80).1
    have hc : ((pdTail : ℚ) : ℝ) = ((log2I (QI.add (QI.pt 1) pdY) 60 80).hi : ℝ)
        + 2 / 41 * (-((log2I pdY 60 80).lo : ℝ)) := by
      unfold pdTail; push_cast <;> ring
    have hT : ((pdTail : ℚ) : ℝ) < 7050097/10^7 := qR_lt_of pdTail_lt (by norm_num)
    linarith
  · by_cases hm6 : m = 6
    · subst hm6
      rw [show ((6:ℕ):ℝ) = (6:ℝ) by norm_num]
      exact le_rfl
    · have key := finite_check m (List.mem_range.mpr (by omega)) hm hm6
      have hs := pdH_sound (m := m) pdY_sound pdY_lo_pos key.1
      have hu : hpd (m:ℝ) ystar ≤ ((pdH m pdY).hi : ℝ) := hs.2
      have hk : ((pdH m pdY).hi : ℝ) < ((pdH 6 pdY).lo : ℝ) := Rat.cast_lt.mpr key.2
      linarith [hW'.1]

theorem pd_constants :
    (16301556/10^7 : ℝ) ≤ (2:ℝ) ^ ωpr ∧ (2:ℝ) ^ ωpr < 16301557/10^7 ∧
    ωpr < 70502/10^5 ∧ (2:ℝ) ^ (70502/10^5 : ℝ) < 16302/10^4 := by
  have hW' : ((pdH 6 pdY).lo : ℝ) ≤ ωpr ∧ ωpr ≤ ((pdH 6 pdY).hi : ℝ) := pdH6_contains
  obtain ⟨hl, hh⟩ := hW'
  have hmono1 : (2:ℝ) ^ ((pdH 6 pdY).lo : ℝ) ≤ (2:ℝ) ^ ωpr :=
    Real.rpow_le_rpow_of_exponent_le (by norm_num) hl
  have hmono2 : (2:ℝ) ^ ωpr ≤ (2:ℝ) ^ ((pdH 6 pdY).hi : ℝ) :=
    Real.rpow_le_rpow_of_exponent_le (by norm_num) hh
  have p1 := (pow2Q_sound (pdH 6 pdY).lo 60 80).1
  have p2 := (pow2Q_sound (pdH 6 pdY).hi 60 80).2
  have p3 := (pow2Q_sound (70502/10^5) 60 80).2
  have hc : ((70502/10^5 : ℚ) : ℝ) = (70502/10^5 : ℝ) := by push_cast <;> ring
  rw [hc] at p3
  have q1 : (16301556/10^7 : ℝ) ≤ ((pow2Q (pdH 6 pdY).lo 60 80).lo : ℝ) :=
    qR_le_of pdC1 (by norm_num)
  have q2 : ((pow2Q (pdH 6 pdY).hi 60 80).hi : ℝ) < 16301557/10^7 :=
    qR_lt_of pdC2 (by norm_num)
  have q3 : ((pdH 6 pdY).hi : ℝ) < 70502/10^5 := qR_lt_of pdC3 (by norm_num)
  have q4 : ((pow2Q (70502/10^5) 60 80).hi : ℝ) < 16302/10^4 :=
    qR_lt_of pdC4 (by norm_num)
  refine ⟨by linarith, by linarith, by linarith, by linarith⟩

end MaxLeafFormal.PDRate