MaxLeafFormal.Enclose

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

import Mathlib

/-!
# Rational enclosures

Computable rational intervals that provably contain `Real.log q`, `Real.logb 2 q`, `Real.exp x`,
`2 ^ x` and `Real.sqrt q` for rational arguments, with interval arithmetic. This is the substitute
for Arb ball arithmetic in the certificate checks; no part of Arb is trusted.

Every function here is defined by structural recursion (fuel) so that `decide +kernel` evaluates
it. Every soundness theorem is unconditional apart from the domain of the real function: when an
internal check fails (for example a poor argument reduction), the functions return a cruder interval
that is still sound.

* `ln (1 - x)`: the series of `Real.abs_log_sub_add_sum_range_le`, after the reduction
  `q = 2^k m` with `m ∈ [1, 2)` and `x = 1 - 1/m`;
* `exp`: the Taylor polynomial with the remainder of `Real.exp_bound`, after halving the argument
  `k` times, then `k` squarings;
* `2^x = exp (x ln 2)`; `logb 2 q = ln q / ln 2`;
* `√q`: a dyadic bracket `ℓ² ≤ q ≤ u²`, checked.

Outward rounding to the denominator `2^prec` keeps the numbers small.
-/

namespace MaxLeafFormal.Enclose

/-- A closed rational interval `[lo, hi]`. -/
structure QI where
  lo : ℚ
  hi : ℚ
deriving DecidableEq, Repr

/-- `I` contains the real `x`. -/
def QI.Contains (I : QI) (x : ℝ) : Prop := (I.lo : ℝ) ≤ x ∧ x ≤ (I.hi : ℝ)

/-! ### Rounding and interval arithmetic -/

/-- Round down to the denominator `2 ^ prec`. -/
def rdn (q : ℚ) (prec : ℕ) : ℚ := ((q * 2 ^ prec).floor : ℚ) / 2 ^ prec

/-- Round up to the denominator `2 ^ prec`. -/
def rup (q : ℚ) (prec : ℕ) : ℚ := ((q * 2 ^ prec).ceil : ℚ) / 2 ^ prec

/-- Outward rounding. -/
def QI.round (I : QI) (prec : ℕ) : QI := ⟨rdn I.lo prec, rup I.hi prec⟩

/-- The point interval. -/
def QI.pt (q : ℚ) : QI := ⟨q, q⟩

def QI.add (I J : QI) : QI := ⟨I.lo + J.lo, I.hi + J.hi⟩

def QI.neg (I : QI) : QI := ⟨-I.hi, -I.lo⟩

def QI.sub (I J : QI) : QI := ⟨I.lo - J.hi, I.hi - J.lo⟩

/-- Scaling by a rational of either sign. -/
def QI.smul (c : ℚ) (I : QI) : QI := if 0 ≤ c then ⟨c * I.lo, c * I.hi⟩ else ⟨c * I.hi, c * I.lo⟩

/-- Product: the extreme corner products. -/
def QI.mul (I J : QI) : QI :=
  ⟨min (min (I.lo * J.lo) (I.lo * J.hi)) (min (I.hi * J.lo) (I.hi * J.hi)),
   max (max (I.lo * J.lo) (I.lo * J.hi)) (max (I.hi * J.lo) (I.hi * J.hi))⟩

/-- Quotient by an interval with positive lower end (soundness needs `0 < J.lo`). -/
def QI.div (I J : QI) : QI :=
  ⟨min (I.lo / J.lo) (I.lo / J.hi), max (I.hi / J.lo) (I.hi / J.hi)⟩

/-! ### Integer helpers with fuel -/

/-- `⌊log₂ n⌋` for `n < 2 ^ fuel` (and `0` for `n ≤ 1`). -/
def ilog2 : ℕ → ℕ → ℕ
  | 0, _ => 0
  | fuel + 1, n => if n < 2 then 0 else ilog2 fuel (n / 2) + 1

/-- Bisection for `⌊√N⌋`: keeps `lo ≤ ⌊√N⌋ < hi`. -/
def isqrtAux : ℕ → ℕ → ℕ → ℕ → ℕ
  | 0, _, lo, _ => lo
  | fuel + 1, N, lo, hi =>
    if hi ≤ lo + 1 then lo
    else
      let m := (lo + hi) / 2
      if m * m ≤ N then isqrtAux fuel N m hi else isqrtAux fuel N lo m

/-- `⌊√N⌋` (a candidate; callers check it). -/
def isqrt (N : ℕ) : ℕ := isqrtAux (2 * ilog2 1024 N + 4) N 0 (N + 1)

/-! ### Logarithm -/

/-- The partial sum `∑_{i<n} x^(i+1)/(i+1)` of `-ln (1 - x)`. -/
def logSeries (x : ℚ) (n : ℕ) : ℚ := ((List.range n).map fun i : ℕ => x ^ (i + 1) / ((i : ℚ) + 1)).sum

/-- The remainder bound `|x|^(n+1) / (1 - |x|)` of `Real.abs_log_sub_add_sum_range_le`. -/
def logErr (x : ℚ) (n : ℕ) : ℚ := |x| ^ (n + 1) / (1 - |x|)

/-- An enclosure of `-ln (1 - x)` for `|x| < 1`. -/
def negLog1mQ (x : ℚ) (n prec : ℕ) : QI :=
  let xl := rdn x prec
  let xh := rup x prec
  if |xl| < 1 ∧ |xh| < 1 then
    ⟨rdn (logSeries xl n - logErr xl n) prec, rup (logSeries xh n + logErr xh n) prec⟩
  else
    -- crude but sound for `x < 1`: `x ≤ -ln (1 - x) ≤ x / (1 - x)`
    ⟨x, if x < 1 then x / (1 - x) else 0⟩

/-- An enclosure of `ln 2 = -ln (1 - 1/2)`. -/
def ln2Q (n prec : ℕ) : QI := negLog1mQ (1 / 2) n prec

/-- The exponent `k` with `q / 2^k ∈ [1, 2)`, for `q > 0` (soundness does not depend on it). -/
def redExp (q : ℚ) : ℤ :=
  let k : ℤ := (ilog2 1024 q.num.natAbs : ℤ) - ilog2 1024 q.den
  if q / 2 ^ k < 1 then k - 1 else k

/-- An enclosure of `ln q`, for `q > 0`. -/
def lnQ (q : ℚ) (n prec : ℕ) : QI :=
  let k := redExp q
  let m := q / 2 ^ k
  let x := 1 - 1 / m
  if 0 ≤ x ∧ x < 1 then
    ((QI.smul (k : ℚ) (ln2Q n prec)).add (negLog1mQ x n prec)).round prec
  else
    -- crude but sound for `q > 0`: `1 - 1/q ≤ ln q ≤ q - 1`
    ⟨1 - 1 / q, q - 1⟩

/-- A positive enclosure of `ln 2` (falling back to `[1/2, 1]`). -/
def ln2Pos (n prec : ℕ) : QI := if 0 < (ln2Q n prec).lo then ln2Q n prec else ⟨1 / 2, 1⟩

/-- An enclosure of `log₂ q = ln q / ln 2`, for `q > 0`. -/
def log2Q (q : ℚ) (n prec : ℕ) : QI := ((lnQ q n prec).div (ln2Pos n prec)).round prec

/-! ### Exponential -/

/-- `∑_{i<n} t^i / i!`. -/
def expSeries (t : ℚ) (n : ℕ) : ℚ :=
  ((List.range n).map fun i => t ^ i / (Nat.factorial i : ℚ)).sum

/-- The remainder bound `|t|^n (n+1) / (n! n)` of `Real.exp_bound` (for `n > 0`). -/
def expErr (t : ℚ) (n : ℕ) : ℚ := |t| ^ n * ((n + 1 : ℚ) / ((Nat.factorial n : ℚ) * n))

/-- Square a nonnegative interval `k` times, rounding outward. -/
def sqIter (prec : ℕ) : ℕ → QI → QI
  | 0, I => I
  | k + 1, I => sqIter prec k ⟨rdn (I.lo * I.lo) prec, rup (I.hi * I.hi) prec⟩

/-- The number of halvings: `|x| / 2^k ≤ 1`. -/
def halvings (x : ℚ) : ℕ := ilog2 1024 (⌈|x|⌉.toNat) + 1

/-- An enclosure of `exp x`. -/
def expQ (x : ℚ) (n prec : ℕ) : QI :=
  let k := halvings x
  let t := x / 2 ^ k
  let N := n + 1
  if |t| ≤ 1 then
    let T := expSeries t N
    let e := expErr t N
    sqIter prec k ⟨max 0 (rdn (T - e) prec), rup (T + e) prec⟩
  else
    -- unreachable in practice; sound: `0 ≤ exp x ≤ 3 ^ ⌈x⌉₊`
    ⟨0, 3 ^ (⌈x⌉.toNat)⟩

/-- An enclosure of `exp` on an interval (monotone image). -/
def expI (I : QI) (n prec : ℕ) : QI := ⟨(expQ I.lo n prec).lo, (expQ I.hi n prec).hi⟩

/-- An enclosure of `2 ^ x = exp (x ln 2)`. -/
def pow2Q (x : ℚ) (n prec : ℕ) : QI := expI ((QI.smul x (ln2Q n prec)).round prec) n prec

/-! ### Square root -/

/-- An enclosure of `√q`, for `q ≥ 0`. -/
def sqrtQ (q : ℚ) (prec : ℕ) : QI :=
  let l : ℚ := (isqrt ⌊q * 4 ^ prec⌋.toNat : ℚ) / 2 ^ prec
  let u : ℚ := l + 1 / 2 ^ prec
  if l * l ≤ q ∧ q ≤ u * u then ⟨l, u⟩ else ⟨0, q + 1⟩

/-! ### Monotone images on intervals -/

/-- `ln` on an interval with positive lower end. -/
def lnI (I : QI) (n prec : ℕ) : QI := ⟨(lnQ I.lo n prec).lo, (lnQ I.hi n prec).hi⟩

/-- `log₂` on an interval with positive lower end. -/
def log2I (I : QI) (n prec : ℕ) : QI := ⟨(log2Q I.lo n prec).lo, (log2Q I.hi n prec).hi⟩

/-- `√` on an interval with nonnegative lower end. -/
def sqrtI (I : QI) (prec : ℕ) : QI := ⟨(sqrtQ I.lo prec).lo, (sqrtQ I.hi prec).hi⟩

/-! ### Soundness -/

section Sound

variable {I J : QI} {x y : ℝ}

lemma enc_mul_mem {a b x y : ℝ} (hx1 : a ≤ x) (hx2 : x ≤ b) :
    min (a * y) (b * y) ≤ x * y ∧ x * y ≤ max (a * y) (b * y) := by
  rcases le_total 0 y with hy | hy
  · exact ⟨(min_le_left _ _).trans (mul_le_mul_of_nonneg_right hx1 hy),
      (mul_le_mul_of_nonneg_right hx2 hy).trans (le_max_right _ _)⟩
  · exact ⟨(min_le_right _ _).trans (mul_le_mul_of_nonpos_right hx2 hy),
      (mul_le_mul_of_nonpos_right hx1 hy).trans (le_max_left _ _)⟩

lemma enc_mul_mem' {c d a y : ℝ} (hy1 : c ≤ y) (hy2 : y ≤ d) :
    min (a * c) (a * d) ≤ a * y ∧ a * y ≤ max (a * c) (a * d) := by
  rcases le_total 0 a with ha | ha
  · exact ⟨(min_le_left _ _).trans (mul_le_mul_of_nonneg_left hy1 ha),
      (mul_le_mul_of_nonneg_left hy2 ha).trans (le_max_right _ _)⟩
  · exact ⟨(min_le_right _ _).trans (mul_le_mul_of_nonpos_left hy2 ha),
      (mul_le_mul_of_nonpos_left hy1 ha).trans (le_max_left _ _)⟩

lemma enc_corner {a b c d x y : ℝ} (hx1 : a ≤ x) (hx2 : x ≤ b) (hy1 : c ≤ y) (hy2 : y ≤ d) :
    min (min (a * c) (a * d)) (min (b * c) (b * d)) ≤ x * y ∧
      x * y ≤ max (max (a * c) (a * d)) (max (b * c) (b * d)) := by
  obtain ⟨h1, h2⟩ := enc_mul_mem (y := y) hx1 hx2
  obtain ⟨ha1, ha2⟩ := enc_mul_mem' (a := a) hy1 hy2
  obtain ⟨hb1, hb2⟩ := enc_mul_mem' (a := b) hy1 hy2
  exact ⟨le_trans (le_min (le_trans (min_le_left _ _) ha1) (le_trans (min_le_right _ _) hb1)) h1,
    le_trans h2 (max_le (le_trans ha2 (le_max_left _ _)) (le_trans hb2 (le_max_right _ _)))⟩

lemma enc_div_lower {a c d y : ℝ} (hc : 0 < c) (hy1 : c ≤ y) (hy2 : y ≤ d) :
    min (a / c) (a / d) ≤ a / y := by
  have hy : 0 < y := lt_of_lt_of_le hc hy1
  have hd : 0 < d := lt_of_lt_of_le hy hy2
  rcases le_total 0 a with ha | ha
  · refine le_trans (min_le_right _ _) ?_
    rw [div_le_div_iff₀ hd hy]
    nlinarith [mul_le_mul_of_nonneg_left hy2 ha]
  · refine le_trans (min_le_left _ _) ?_
    rw [div_le_div_iff₀ hc hy]
    nlinarith [mul_le_mul_of_nonpos_left hy1 ha]

lemma enc_div_upper {a c d y : ℝ} (hc : 0 < c) (hy1 : c ≤ y) (hy2 : y ≤ d) :
    a / y ≤ max (a / c) (a / d) := by
  have hy : 0 < y := lt_of_lt_of_le hc hy1
  have hd : 0 < d := lt_of_lt_of_le hy hy2
  rcases le_total 0 a with ha | ha
  · refine le_trans ?_ (le_max_left _ _)
    rw [div_le_div_iff₀ hy hc]
    nlinarith [mul_le_mul_of_nonneg_left hy1 ha]
  · refine le_trans ?_ (le_max_right _ _)
    rw [div_le_div_iff₀ hy hd]
    nlinarith [mul_le_mul_of_nonpos_left hy2 ha]

lemma enc_sqrt_bracket {q l u : ℚ} (hl : 0 ≤ l) (hu : 0 ≤ u) (h : l * l ≤ q ∧ q ≤ u * u) :
    (l : ℝ) ≤ Real.sqrt q ∧ Real.sqrt q ≤ (u : ℝ) := by
  obtain ⟨h1, h2⟩ := h
  constructor
  · have h1' : ((l * l : ℚ) : ℝ) ≤ (q : ℝ) := by exact_mod_cast h1
    push_cast at h1'
    have := Real.sqrt_le_sqrt h1'
    rwa [Real.sqrt_mul_self (show (0:ℝ) ≤ (l:ℝ) by exact_mod_cast hl)] at this
  · have h2' : (q : ℝ) ≤ ((u * u : ℚ) : ℝ) := by exact_mod_cast h2
    push_cast at h2'
    have := Real.sqrt_le_sqrt h2'
    rwa [Real.sqrt_mul_self (show (0:ℝ) ≤ (u:ℝ) by exact_mod_cast hu)] at this

theorem rdn_le (q : ℚ) (prec : ℕ) : rdn q prec ≤ q := by
  unfold rdn
  have h2 : (0:ℚ) < 2 ^ prec := by positivity
  rw [div_le_iff₀ h2]
  exact Int.floor_le (q * 2 ^ prec)

theorem le_rup (q : ℚ) (prec : ℕ) : q ≤ rup q prec := by
  unfold rup
  have h2 : (0:ℚ) < 2 ^ prec := by positivity
  rw [le_div_iff₀ h2]
  exact Rat.le_ceil

theorem QI.round_sound (h : I.Contains x) (prec : ℕ) : (I.round prec).Contains x := by
  simp only [QI.Contains] at h ⊢
  obtain ⟨h1, h2⟩ := h
  simp only [QI.round]
  have e1 : ((rdn I.lo prec : ℚ) : ℝ) ≤ (I.lo : ℝ) := by exact_mod_cast rdn_le I.lo prec
  have e2 : (I.hi : ℝ) ≤ ((rup I.hi prec : ℚ) : ℝ) := by exact_mod_cast le_rup I.hi prec
  exact ⟨le_trans e1 h1, le_trans h2 e2⟩

theorem QI.pt_sound (q : ℚ) : (QI.pt q).Contains (q : ℝ) := by
  simp [QI.Contains, QI.pt]

theorem QI.add_sound (hI : I.Contains x) (hJ : J.Contains y) : (I.add J).Contains (x + y) := by
  simp only [QI.Contains] at hI hJ ⊢
  obtain ⟨h1, h2⟩ := hI
  obtain ⟨k1, k2⟩ := hJ
  simp only [QI.add]
  push_cast
  exact ⟨add_le_add h1 k1, add_le_add h2 k2⟩

theorem QI.neg_sound (hI : I.Contains x) : I.neg.Contains (-x) := by
  simp only [QI.Contains] at hI ⊢
  obtain ⟨h1, h2⟩ := hI
  simp only [QI.neg]
  push_cast
  exact ⟨neg_le_neg h2, neg_le_neg h1⟩

theorem QI.sub_sound (hI : I.Contains x) (hJ : J.Contains y) : (I.sub J).Contains (x - y) := by
  simp only [QI.Contains] at hI hJ ⊢
  obtain ⟨h1, h2⟩ := hI
  obtain ⟨k1, k2⟩ := hJ
  simp only [QI.sub]
  push_cast
  exact ⟨sub_le_sub h1 k2, sub_le_sub h2 k1⟩

theorem QI.smul_sound (c : ℚ) (hI : I.Contains x) : (I.smul c).Contains ((c : ℝ) * x) := by
  simp only [QI.Contains] at hI ⊢
  obtain ⟨h1, h2⟩ := hI
  by_cases h : 0 ≤ c
  · have hc : (0:ℝ) ≤ c := by exact_mod_cast h
    simp only [QI.smul, ite_eq_left h]
    push_cast
    exact ⟨mul_le_mul_of_nonneg_left h1 hc, mul_le_mul_of_nonneg_left h2 hc⟩
  · have h' : c ≤ 0 := (not_le.mp h).le
    have hc : (c:ℝ) ≤ 0 := by exact_mod_cast h'
    simp only [QI.smul, ite_eq_right h]
    push_cast
    exact ⟨mul_le_mul_of_nonpos_left h2 hc, mul_le_mul_of_nonpos_left h1 hc⟩

theorem QI.mul_sound (hI : I.Contains x) (hJ : J.Contains y) : (I.mul J).Contains (x * y) := by
  simp only [QI.Contains] at hI hJ ⊢
  obtain ⟨h1, h2⟩ := hI
  obtain ⟨k1, k2⟩ := hJ
  simp only [QI.mul]
  push_cast
  exact enc_corner h1 h2 k1 k2

theorem QI.div_sound (hI : I.Contains x) (hJ : J.Contains y) (hJ0 : 0 < J.lo) :
    (I.div J).Contains (x / y) := by
  simp only [QI.Contains] at hI hJ ⊢
  obtain ⟨h1, h2⟩ := hI
  obtain ⟨k1, k2⟩ := hJ
  have hc : (0:ℝ) < (J.lo:ℝ) := by exact_mod_cast hJ0
  have hy : 0 < y := lt_of_lt_of_le hc k1
  simp only [QI.div]
  push_cast
  constructor
  · exact le_trans (enc_div_lower (a := (I.lo:ℝ)) hc k1 k2) (div_le_div_of_nonneg_right h1 hy.le)
  · exact le_trans (div_le_div_of_nonneg_right h2 hy.le) (enc_div_upper (a := (I.hi:ℝ)) hc k1 k2)

theorem logSeries_succ (x : ℚ) (n : ℕ) :
    logSeries x (n + 1) = logSeries x n + x ^ (n + 1) / ((n : ℚ) + 1) := by
  unfold logSeries
  rw [List.range_succ, List.map_append, List.sum_append]
  simp

theorem logSeries_cast (x : ℚ) (n : ℕ) :
    ((logSeries x n : ℚ) : ℝ) = ∑ i ∈ Finset.range n, (x : ℝ) ^ (i + 1) / ((i : ℝ) + 1) := by
  induction n with
  | zero => simp [logSeries]
  | succ n ih =>
    rw [logSeries_succ, Rat.cast_add, ih, Finset.sum_range_succ]
    push_cast <;> ring

theorem logErr_cast (x : ℚ) (n : ℕ) :
    ((logErr x n : ℚ) : ℝ) = |(x : ℝ)| ^ (n + 1) / (1 - |(x : ℝ)|) := by
  unfold logErr
  push_cast <;> ring

theorem expSeries_succ (t : ℚ) (n : ℕ) :
    expSeries t (n + 1) = expSeries t n + t ^ n / (Nat.factorial n : ℚ) := by
  unfold expSeries
  rw [List.range_succ, List.map_append, List.sum_append]
  simp

theorem expSeries_cast (t : ℚ) (n : ℕ) :
    ((expSeries t n : ℚ) : ℝ) = ∑ i ∈ Finset.range n, (t : ℝ) ^ i / ((Nat.factorial i : ℕ) : ℝ) := by
  induction n with
  | zero => simp [expSeries]
  | succ n ih =>
    rw [expSeries_succ, Rat.cast_add, ih, Finset.sum_range_succ]
    push_cast <;> ring

theorem expErr_cast (t : ℚ) (N : ℕ) :
    ((expErr t N : ℚ) : ℝ) =
      |(t : ℝ)| ^ N * (((N.succ : ℕ) : ℝ) / (((Nat.factorial N : ℕ) : ℝ) * (N : ℝ))) := by
  unfold expErr
  push_cast <;> ring

theorem negLog_lower {a X : ℝ} (n : ℕ) (ha : |a| < 1) (hax : a ≤ X) (hX : X < 1) :
    (∑ i ∈ Finset.range n, a ^ (i + 1) / ((i : ℝ) + 1)) - |a| ^ (n + 1) / (1 - |a|)
      ≤ -Real.log (1 - X) := by
  have h1 := (abs_le.mp (Real.abs_log_sub_add_sum_range_le ha n)).2
  have hlog : Real.log (1 - X) ≤ Real.log (1 - a) :=
    Real.log_le_log (by linarith) (by linarith)
  linarith

theorem negLog_upper {b X : ℝ} (n : ℕ) (hb : |b| < 1) (hXb : X ≤ b) :
    -Real.log (1 - X) ≤
      (∑ i ∈ Finset.range n, b ^ (i + 1) / ((i : ℝ) + 1)) + |b| ^ (n + 1) / (1 - |b|) := by
  have h1 := (abs_le.mp (Real.abs_log_sub_add_sum_range_le hb n)).1
  have hb1 : b < 1 := (abs_lt.mp hb).2
  have hlog : Real.log (1 - b) ≤ Real.log (1 - X) :=
    Real.log_le_log (by linarith) (by linarith)
  linarith

theorem negLog1mQ_sound {x : ℚ} (hx : x < 1) (n prec : ℕ) :
    (negLog1mQ x n prec).Contains (-Real.log (1 - x)) := by
  unfold negLog1mQ
  dsimp only
  by_cases h : |rdn x prec| < 1 ∧ |rup x prec| < 1
  · rw [ite_eq_left h]
    obtain ⟨hl, hh⟩ := h
    have hA : |((rdn x prec : ℚ) : ℝ)| < 1 := by
      rw [← Rat.cast_abs]; exact_mod_cast hl
    have hB : |((rup x prec : ℚ) : ℝ)| < 1 := by
      rw [← Rat.cast_abs]; exact_mod_cast hh
    have ha1 : ((rdn x prec : ℚ) : ℝ) ≤ ((x : ℚ) : ℝ) := Rat.cast_le.mpr (rdn_le x prec)
    have hb1 : ((x : ℚ) : ℝ) ≤ ((rup x prec : ℚ) : ℝ) := Rat.cast_le.mpr (le_rup x prec)
    have hX : ((x : ℚ) : ℝ) < 1 := by exact_mod_cast hx
    have hL := negLog_lower (a := ((rdn x prec : ℚ) : ℝ)) (X := ((x : ℚ) : ℝ)) n hA ha1 hX
    have hU := negLog_upper (b := ((rup x prec : ℚ) : ℝ)) (X := ((x : ℚ) : ℝ)) n hB hb1
    have c1 : ((logSeries (rdn x prec) n - logErr (rdn x prec) n : ℚ) : ℝ)
        ≤ -Real.log (1 - ((x : ℚ) : ℝ)) := by
      rw [Rat.cast_sub, logSeries_cast, logErr_cast]; exact hL
    have c2 : -Real.log (1 - ((x : ℚ) : ℝ))
        ≤ ((logSeries (rup x prec) n + logErr (rup x prec) n : ℚ) : ℝ) := by
      rw [Rat.cast_add, logSeries_cast, logErr_cast]; exact hU
    unfold QI.Contains
    exact ⟨le_trans (Rat.cast_le.mpr (rdn_le _ _)) c1, le_trans c2 (Rat.cast_le.mpr (le_rup _ _))⟩
  · rw [ite_eq_right h, ite_eq_left hx]
    have hX : ((x : ℚ) : ℝ) < 1 := by exact_mod_cast hx
    have hp : (0 : ℝ) < 1 - ((x : ℚ) : ℝ) := by linarith
    have hne : (1 - ((x : ℚ) : ℝ)) ≠ 0 := hp.ne'
    have l1 := Real.log_le_sub_one_of_pos hp
    have l2 := Real.one_sub_inv_le_log_of_pos hp
    have hi := mul_inv_cancel₀ hne
    have e : ((x : ℚ) : ℝ) / (1 - ((x : ℚ) : ℝ)) = (1 - ((x : ℚ) : ℝ))⁻¹ - 1 := by
      rw [div_eq_mul_inv]
      linear_combination (-1 : ℝ) * hi
    unfold QI.Contains
    refine ⟨?_, ?_⟩
    · show ((x : ℚ) : ℝ) ≤ -Real.log (1 - ((x : ℚ) : ℝ))
      linarith
    · show -Real.log (1 - ((x : ℚ) : ℝ)) ≤ ((x / (1 - x) : ℚ) : ℝ)
      push_cast
      rw [e]
      linarith

theorem ln2Q_sound (n prec : ℕ) : (ln2Q n prec).Contains (Real.log 2) := by
  have h := negLog1mQ_sound (x := 1 / 2) (by norm_num) n prec
  have e : -Real.log (1 - ((1 / 2 : ℚ) : ℝ)) = Real.log 2 := by
    rw [show (1 - ((1 / 2 : ℚ) : ℝ)) = (2 : ℝ)⁻¹ by norm_num, Real.log_inv, neg_neg]
  unfold ln2Q
  rw [← e]
  exact h

theorem log_split (q : ℚ) (hq : 0 < q) (k : ℤ) :
    Real.log (q : ℝ) = (k : ℝ) * Real.log 2 + Real.log ((q : ℝ) / 2 ^ k) := by
  have hq' : (0 : ℝ) < q := Rat.cast_pos.mpr hq
  have h2 : (0 : ℝ) < 2 ^ k := by positivity
  rw [Real.log_div hq'.ne' h2.ne', Real.log_zpow]
  ring

theorem neg_log_one_sub_inv (m : ℝ) : -Real.log (1 - (1 - 1 / m)) = Real.log m := by
  rw [show (1 - (1 - 1 / m)) = m⁻¹ by ring, Real.log_inv, neg_neg]

theorem lnQ_main (q : ℚ) (hq : 0 < q) (k : ℤ) (n prec : ℕ)
    (h : 0 ≤ 1 - 1 / (q / 2 ^ k) ∧ 1 - 1 / (q / 2 ^ k) < 1) :
    (((QI.smul (k : ℚ) (ln2Q n prec)).add (negLog1mQ (1 - 1 / (q / 2 ^ k)) n prec)).round
      prec).Contains (Real.log q) := by
  have hlog := log_split q hq k
  have hneg := negLog1mQ_sound h.2 n prec
  have hc : ((1 - 1 / (q / 2 ^ k) : ℚ) : ℝ) = 1 - 1 / ((q : ℝ) / 2 ^ k) := by
    push_cast <;> ring
  rw [hc, neg_log_one_sub_inv] at hneg
  have hsm := QI.smul_sound (k : ℚ) (ln2Q_sound n prec)
  have hsum := QI.add_sound hsm hneg
  have e : ((k : ℚ) : ℝ) * Real.log 2 + Real.log ((q : ℝ) / 2 ^ k) = Real.log q := by
    rw [hlog, Rat.cast_intCast]
  rw [← e]
  exact QI.round_sound hsum prec

theorem lnQ_sound {q : ℚ} (hq : 0 < q) (n prec : ℕ) : (lnQ q n prec).Contains (Real.log q) := by
  unfold lnQ
  dsimp only
  split_ifs with h
  · exact lnQ_main q hq _ n prec h
  · have hq' : (0 : ℝ) < q := Rat.cast_pos.mpr hq
    unfold QI.Contains
    refine ⟨?_, ?_⟩
    · show ((1 - 1 / q : ℚ) : ℝ) ≤ Real.log q
      push_cast
      have := Real.one_sub_inv_le_log_of_pos hq'
      rw [one_div]
      linarith
    · show Real.log q ≤ ((q - 1 : ℚ) : ℝ)
      push_cast
      linarith [Real.log_le_sub_one_of_pos hq']

theorem ln2Pos_sound' (n prec : ℕ) :
    (ln2Pos n prec).Contains (Real.log 2) ∧ 0 < (ln2Pos n prec).lo := by
  unfold ln2Pos
  split_ifs with h
  · exact ⟨ln2Q_sound n prec, h⟩
  · refine ⟨?_, by show (0 : ℚ) < 1 / 2; norm_num⟩
    unfold QI.Contains
    constructor
    · show ((1 / 2 : ℚ) : ℝ) ≤ Real.log 2
      have := Real.log_two_gt_d9
      norm_num at this ⊢
      linarith
    · show Real.log 2 ≤ ((1 : ℚ) : ℝ)
      have := Real.log_two_lt_d9
      norm_num at this ⊢
      linarith

theorem log2Q_sound {q : ℚ} (hq : 0 < q) (n prec : ℕ) :
    (log2Q q n prec).Contains (Real.logb 2 q) := by
  have hD := QI.div_sound (lnQ_sound hq n prec) (ln2Pos_sound' n prec).1 (ln2Pos_sound' n prec).2
  unfold log2Q
  show (((lnQ q n prec).div (ln2Pos n prec)).round prec).Contains (Real.log q / Real.log 2)
  exact QI.round_sound hD prec

theorem rdn_nonneg {q : ℚ} (hq : 0 ≤ q) (prec : ℕ) : 0 ≤ rdn q prec := by
  unfold rdn
  refine div_nonneg ?_ (by positivity)
  exact_mod_cast Int.floor_nonneg.mpr (mul_nonneg hq (by positivity))

theorem sqIter_sound (prec k : ℕ) :
    ∀ (I : QI) (v : ℝ), I.Contains v → 0 ≤ I.lo → (sqIter prec k I).Contains (v ^ (2 ^ k)) := by
  induction k with
  | zero =>
    intro I v h _
    simpa [sqIter] using h
  | succ k ih =>
    intro I v h h0
    simp only [QI.Contains] at h
    obtain ⟨h1, h2⟩ := h
    have hl0 : (0 : ℝ) ≤ I.lo := by exact_mod_cast h0
    have hv0 : (0 : ℝ) ≤ v := hl0.trans h1
    have hsq1 : (I.lo : ℝ) * I.lo ≤ v * v := mul_le_mul h1 h1 hl0 hv0
    have hsq2 : v * v ≤ (I.hi : ℝ) * I.hi := mul_le_mul h2 h2 hv0 (hv0.trans h2)
    have hA : ((rdn (I.lo * I.lo) prec : ℚ) : ℝ) ≤ v ^ 2 := by
      have : ((rdn (I.lo * I.lo) prec : ℚ) : ℝ) ≤ ((I.lo * I.lo : ℚ) : ℝ) :=
        Rat.cast_le.mpr (rdn_le (I.lo * I.lo) prec)
      push_cast at this
      nlinarith
    have hB : v ^ 2 ≤ ((rup (I.hi * I.hi) prec : ℚ) : ℝ) := by
      have : ((I.hi * I.hi : ℚ) : ℝ) ≤ ((rup (I.hi * I.hi) prec : ℚ) : ℝ) :=
        Rat.cast_le.mpr (le_rup (I.hi * I.hi) prec)
      push_cast at this
      nlinarith
    have hJ0 : 0 ≤ rdn (I.lo * I.lo) prec := rdn_nonneg (mul_nonneg h0 h0) prec
    have key := ih ⟨rdn (I.lo * I.lo) prec, rup (I.hi * I.hi) prec⟩ (v ^ 2)
      (by unfold QI.Contains; exact ⟨hA, hB⟩) hJ0
    have e : (v ^ 2) ^ (2 ^ k) = v ^ (2 ^ (k + 1)) := by
      rw [← pow_mul, pow_succ, mul_comm]
    rw [e] at key
    simp only [sqIter]
    exact key

theorem max_zero_cast_le {r : ℚ} {v : ℝ} (hr : (r : ℝ) ≤ v) (hv : 0 ≤ v) :
    ((max 0 r : ℚ) : ℝ) ≤ v := by
  rcases le_total 0 r with h | h
  · rw [max_eq_right h]; exact hr
  · rw [max_eq_left h]; simpa using hv

theorem abs_cast_le {t c : ℚ} (h : |t| ≤ c) : |(t : ℝ)| ≤ (c : ℝ) := by
  have : ((|t| : ℚ) : ℝ) ≤ ((c : ℚ) : ℝ) := Rat.cast_le.mpr h
  push_cast at this
  exact this

theorem exp_core (x t : ℚ) (k n prec : ℕ) (hxt : x = t * 2 ^ k) (ht : |t| ≤ 1) :
    (sqIter prec k ⟨max 0 (rdn (expSeries t (n + 1) - expErr t (n + 1)) prec),
      rup (expSeries t (n + 1) + expErr t (n + 1)) prec⟩).Contains (Real.exp x) := by
  have hR : |(t : ℝ)| ≤ 1 := by
    have := abs_cast_le ht
    rwa [Rat.cast_one] at this
  have hS := Real.exp_bound hR (n := n + 1) (Nat.succ_pos n)
  have hb := abs_le.mp hS
  have hT := expSeries_cast t (n + 1)
  have hE := expErr_cast t (n + 1)
  have hlo : ((expSeries t (n + 1) - expErr t (n + 1) : ℚ) : ℝ) ≤ Real.exp t := by
    rw [Rat.cast_sub, hT, hE]; linarith [hb.1]
  have hhi : Real.exp t ≤ ((expSeries t (n + 1) + expErr t (n + 1) : ℚ) : ℝ) := by
    rw [Rat.cast_add, hT, hE]; linarith [hb.2]
  have hI : (⟨max 0 (rdn (expSeries t (n + 1) - expErr t (n + 1)) prec),
      rup (expSeries t (n + 1) + expErr t (n + 1)) prec⟩ : QI).Contains (Real.exp t) := by
    unfold QI.Contains
    exact ⟨max_zero_cast_le (le_trans (Rat.cast_le.mpr (rdn_le _ _)) hlo) (Real.exp_pos _).le,
      le_trans hhi (Rat.cast_le.mpr (le_rup _ _))⟩
  have key := sqIter_sound prec k _ (Real.exp t) hI (le_max_left _ _)
  have hexp : Real.exp (x : ℝ) = (Real.exp (t : ℝ)) ^ (2 ^ k) := by
    rw [← Real.exp_nat_mul]
    congr 1
    rw [hxt]
    push_cast <;> ring
  rw [hexp]
  exact key

theorem expQ_sound (x : ℚ) (n prec : ℕ) : (expQ x n prec).Contains (Real.exp x) := by
  unfold expQ
  dsimp only
  split_ifs with h
  · exact exp_core x (x / 2 ^ halvings x) (halvings x) n prec (by field_simp) h
  · unfold QI.Contains
    refine ⟨?_, ?_⟩
    · show ((0 : ℚ) : ℝ) ≤ Real.exp x
      rw [Rat.cast_zero]
      exact (Real.exp_pos _).le
    · show Real.exp x ≤ ((3 ^ (⌈x⌉.toNat) : ℚ) : ℝ)
      push_cast
      have h1 : x ≤ ((⌈x⌉ : ℤ) : ℚ) := Int.le_ceil x
      have h2 : (⌈x⌉ : ℤ) ≤ ((⌈x⌉.toNat : ℕ) : ℤ) := Int.self_le_toNat _
      have h3 : (x : ℝ) ≤ ((⌈x⌉ : ℤ) : ℝ) := by exact_mod_cast h1
      have h4 : ((⌈x⌉ : ℤ) : ℝ) ≤ ((⌈x⌉.toNat : ℕ) : ℝ) := by exact_mod_cast h2
      have hx : (x : ℝ) ≤ ((⌈x⌉.toNat : ℕ) : ℝ) := h3.trans h4
      have he : Real.exp 1 ≤ 3 := by
        have := Real.exp_one_lt_d9
        norm_num at this
        linarith
      calc Real.exp x ≤ Real.exp ((⌈x⌉.toNat : ℕ) : ℝ) := Real.exp_le_exp.mpr hx
        _ = Real.exp 1 ^ (⌈x⌉.toNat) := by rw [← Real.exp_nat_mul, mul_one]
        _ ≤ 3 ^ (⌈x⌉.toNat) := pow_le_pow_left₀ (Real.exp_pos 1).le he _

theorem expI_sound (hI : I.Contains x) (n prec : ℕ) : (expI I n prec).Contains (Real.exp x) := by
  simp only [QI.Contains] at hI ⊢
  obtain ⟨hl, hh⟩ := hI
  have h1 := expQ_sound I.lo n prec
  have h2 := expQ_sound I.hi n prec
  simp only [QI.Contains] at h1 h2
  unfold expI
  exact ⟨le_trans h1.1 (Real.exp_le_exp.mpr hl), le_trans (Real.exp_le_exp.mpr hh) h2.2⟩

theorem pow2Q_sound (x : ℚ) (n prec : ℕ) : (pow2Q x n prec).Contains ((2:ℝ) ^ (x:ℝ)) := by
  have hs := QI.smul_sound x (ln2Q_sound n prec)
  have hr := QI.round_sound hs prec
  have he := expI_sound hr n prec
  unfold pow2Q
  rw [Real.rpow_def_of_pos (by norm_num : (0:ℝ) < 2), mul_comm (Real.log 2)]
  exact he

theorem sqrtQ_sound {q : ℚ} (hq : 0 ≤ q) (prec : ℕ) : (sqrtQ q prec).Contains (Real.sqrt q) := by
  simp only [sqrtQ]
  split_ifs with h
  · simp only [QI.Contains]
    exact enc_sqrt_bracket (by positivity) (by positivity) h
  · simp only [QI.Contains]
    push_cast
    have hqr : (0:ℝ) ≤ (q:ℝ) := by exact_mod_cast hq
    refine ⟨Real.sqrt_nonneg _, ?_⟩
    have := Real.sqrt_le_sqrt (show (q:ℝ) ≤ ((q:ℝ) + 1)^2 by nlinarith)
    rwa [Real.sqrt_sq (show (0:ℝ) ≤ (q:ℝ) + 1 by linarith)] at this

theorem lnI_sound (hI : I.Contains x) (h0 : 0 < I.lo) (n prec : ℕ) :
    (lnI I n prec).Contains (Real.log x) := by
  simp only [QI.Contains] at hI
  obtain ⟨hl, hh⟩ := hI
  have hlh : I.lo ≤ I.hi := by exact_mod_cast hl.trans hh
  have hlo : (0 : ℝ) < I.lo := by exact_mod_cast h0
  have hx : (0 : ℝ) < x := lt_of_lt_of_le hlo hl
  have h1 := lnQ_sound h0 n prec
  have h2 := lnQ_sound (lt_of_lt_of_le h0 hlh) n prec
  simp only [QI.Contains] at h1 h2
  unfold lnI
  simp only [QI.Contains]
  exact ⟨le_trans h1.1 (Real.log_le_log hlo hl), le_trans (Real.log_le_log hx hh) h2.2⟩

theorem log2I_sound (hI : I.Contains x) (h0 : 0 < I.lo) (n prec : ℕ) :
    (log2I I n prec).Contains (Real.logb 2 x) := by
  simp only [QI.Contains] at hI
  obtain ⟨hl, hh⟩ := hI
  have hlh : I.lo ≤ I.hi := by exact_mod_cast hl.trans hh
  have hlo : (0 : ℝ) < I.lo := by exact_mod_cast h0
  have hx : (0 : ℝ) < x := lt_of_lt_of_le hlo hl
  have h1 := log2Q_sound h0 n prec
  have h2 := log2Q_sound (lt_of_lt_of_le h0 hlh) n prec
  simp only [QI.Contains] at h1 h2
  unfold log2I
  simp only [QI.Contains]
  exact ⟨le_trans h1.1 (Real.logb_le_logb_of_le (by norm_num) hlo hl),
    le_trans (Real.logb_le_logb_of_le (by norm_num) hx hh) h2.2⟩

theorem sqrtI_sound (hI : I.Contains x) (h0 : 0 ≤ I.lo) (prec : ℕ) :
    (sqrtI I prec).Contains (Real.sqrt x) := by
  simp only [QI.Contains] at hI
  obtain ⟨h1, h2⟩ := hI
  have hlr : (I.lo : ℝ) ≤ I.hi := h1.trans h2
  have hhi : (0:ℚ) ≤ I.hi := le_trans h0 (by exact_mod_cast hlr)
  have hlo := sqrtQ_sound h0 prec
  have hhi' := sqrtQ_sound hhi prec
  obtain ⟨a1, a2⟩ := hlo
  obtain ⟨b1, b2⟩ := hhi'
  have c1 : ((sqrtQ I.lo prec).lo : ℝ) ≤ Real.sqrt x := le_trans a1 (Real.sqrt_le_sqrt h1)
  have c2 : Real.sqrt x ≤ ((sqrtQ I.hi prec).hi : ℝ) := le_trans (Real.sqrt_le_sqrt h2) b2
  unfold sqrtI QI.Contains
  exact ⟨c1, c2⟩

end Sound

end MaxLeafFormal.Enclose