IndisputableMonolith.Numerics.Interval.Trig

   1import IndisputableMonolith.Numerics.Interval.Basic
   2import IndisputableMonolith.Numerics.Interval.PiBounds
   3import Mathlib.Analysis.SpecialFunctions.Trigonometric.Arctan
   4import Mathlib.Analysis.SpecialFunctions.Trigonometric.ArctanDeriv
   5import Mathlib.Analysis.SpecialFunctions.Complex.Arctan
   6import Mathlib.Analysis.SpecificLimits.Normed
   7import Mathlib.Analysis.Calculus.MeanValue
   8
   9/-!
  10# Interval Arithmetic for Trigonometric Functions
  11
  12This module provides rigorous interval bounds for arctan and tan.
  13
  14## Constructive proofs
  15
  16All arctan bounds are proved constructively using derivative-comparison
  17monotonicity (`monotoneOn_of_deriv_nonneg`):
  18
  19- For the **upper bound**: `arctan(x) ≤ x − x³/3 + x⁵/5`, because
  20  `(1−t²+t⁴) ≥ 1/(1+t²)` (equivalently `1+t⁶ ≥ 1`).
  21
  22- For the **lower bound**: `x − x³/3 + x⁵/5 − x⁷/7 ≤ arctan(x)`, because
  23  `(1−t²+t⁴−t⁶) ≤ 1/(1+t²)` (equivalently `1−t⁸ ≤ 1`).
  24
  25Both polynomials are antiderivatives of their respective bounding functions,
  26and the comparison follows from `monotoneOn_of_deriv_nonneg` on `[0,∞)`.
  27
  28The arctan(2) identity `arctan(2) = π/4 + arctan(1/3)` is proved using
  29`Real.arctan_add` from Mathlib.
  30-/
  31
  32namespace IndisputableMonolith.Numerics
  33
  34open Interval
  35open Real
  36open Filter
  37open scoped Topology Finset
  38
  39/-! ## §1. Taylor-polynomial bounds for arctan via derivative comparison -/
  40
  41/-- Upper bounding polynomial: g(x) = x − x³/3 + x⁵/5 -/
  42private noncomputable def g_upper (x : ℝ) : ℝ := x - x ^ 3 / 3 + x ^ 5 / 5
  43
  44private theorem g_upper_continuous : Continuous g_upper := by unfold g_upper; fun_prop
  45private theorem g_upper_differentiable : Differentiable ℝ g_upper := by unfold g_upper; fun_prop
  46
  47private theorem g_upper_deriv (t : ℝ) :
  48    HasDerivAt g_upper (1 - t ^ 2 + t ^ 4) t := by
  49  unfold g_upper
  50  have := ((hasDerivAt_id t).sub ((hasDerivAt_pow 3 t).div_const 3)).add
  51    ((hasDerivAt_pow 5 t).div_const 5)
  52  convert this using 1; ring
  53
  54/-- Key inequality: `1/(1+t²) ≤ 1 − t² + t⁴` for all t.
  55    Proof: `(1−t²+t⁴)(1+t²) = 1+t⁶ ≥ 1`. -/
  56private theorem inv_one_add_sq_le_upper (t : ℝ) :
  57    1 / (1 + t ^ 2) ≤ 1 - t ^ 2 + t ^ 4 := by
  58  rw [div_le_iff₀ (by positivity : 0 < 1 + t ^ 2)]
  59  have : (1 - t ^ 2 + t ^ 4) * (1 + t ^ 2) = 1 + t ^ 6 := by ring
  60  rw [this]; linarith [sq_nonneg (t ^ 3)]
  61
  62/-- `arctan(x) ≤ x − x³/3 + x⁵/5` for x ≥ 0. -/
  63theorem arctan_le_upper_poly (x : ℝ) (hx : 0 ≤ x) : arctan x ≤ g_upper x := by
  64  suffices h : 0 ≤ g_upper x - arctan x by linarith
  65  have hkey : MonotoneOn (fun t => g_upper t - arctan t) (Set.Ici 0) :=
  66    monotoneOn_of_deriv_nonneg (convex_Ici 0)
  67      ((g_upper_continuous.sub continuous_arctan).continuousOn)
  68      (fun t _ => ((g_upper_differentiable t).sub
  69        (hasDerivAt_arctan t).differentiableAt).differentiableWithinAt)
  70      (fun t ht => by
  71        simp only [Set.nonempty_Iio, interior_Ici'] at ht
  72        have hd : HasDerivAt (fun s => g_upper s - arctan s)
  73          ((1 - t^2 + t^4) - 1/(1+t^2)) t :=
  74          (g_upper_deriv t).sub (hasDerivAt_arctan t)
  75        rw [hd.deriv]
  76        linarith [inv_one_add_sq_le_upper t])
  77  linarith [hkey (Set.mem_Ici.mpr (le_refl 0)) (Set.mem_Ici.mpr hx) hx,
  78            show g_upper 0 - arctan 0 = 0 by simp [g_upper, arctan_zero]]
  79
  80/-- Lower bounding polynomial: h(x) = x − x³/3 + x⁵/5 − x⁷/7 -/
  81private noncomputable def h_lower (x : ℝ) : ℝ := x - x ^ 3 / 3 + x ^ 5 / 5 - x ^ 7 / 7
  82
  83private theorem h_lower_continuous : Continuous h_lower := by unfold h_lower; fun_prop
  84private theorem h_lower_differentiable : Differentiable ℝ h_lower := by unfold h_lower; fun_prop
  85
  86private theorem h_lower_deriv (t : ℝ) :
  87    HasDerivAt h_lower (1 - t ^ 2 + t ^ 4 - t ^ 6) t := by
  88  unfold h_lower
  89  have := (((hasDerivAt_id t).sub ((hasDerivAt_pow 3 t).div_const 3)).add
  90    ((hasDerivAt_pow 5 t).div_const 5)).sub ((hasDerivAt_pow 7 t).div_const 7)
  91  convert this using 1; ring
  92
  93/-- Key inequality: `1 − t² + t⁴ − t⁶ ≤ 1/(1+t²)` for all t.
  94    Proof: `(1−t²+t⁴−t⁶)(1+t²) = 1−t⁸ ≤ 1`. -/
  95private theorem lower_le_inv_one_add_sq (t : ℝ) :
  96    1 - t ^ 2 + t ^ 4 - t ^ 6 ≤ 1 / (1 + t ^ 2) := by
  97  rw [le_div_iff₀ (by positivity : 0 < 1 + t ^ 2)]
  98  have : (1 - t ^ 2 + t ^ 4 - t ^ 6) * (1 + t ^ 2) = 1 - t ^ 8 := by ring
  99  rw [this]; linarith [sq_nonneg (t ^ 4)]
 100
 101/-- `x − x³/3 + x⁵/5 − x⁷/7 ≤ arctan(x)` for x ≥ 0. -/
 102theorem lower_poly_le_arctan (x : ℝ) (hx : 0 ≤ x) : h_lower x ≤ arctan x := by
 103  suffices h : 0 ≤ arctan x - h_lower x by linarith
 104  have hkey : MonotoneOn (fun t => arctan t - h_lower t) (Set.Ici 0) :=
 105    monotoneOn_of_deriv_nonneg (convex_Ici 0)
 106      ((continuous_arctan.sub h_lower_continuous).continuousOn)
 107      (fun t _ => ((hasDerivAt_arctan t).differentiableAt.sub
 108        (h_lower_differentiable t)).differentiableWithinAt)
 109      (fun t ht => by
 110        simp only [Set.nonempty_Iio, interior_Ici'] at ht
 111        have hd : HasDerivAt (fun s => arctan s - h_lower s)
 112          (1/(1+t^2) - (1 - t^2 + t^4 - t^6)) t :=
 113          (hasDerivAt_arctan t).sub (h_lower_deriv t)
 114        rw [hd.deriv]
 115        linarith [lower_le_inv_one_add_sq t])
 116  linarith [hkey (Set.mem_Ici.mpr (le_refl 0)) (Set.mem_Ici.mpr hx) hx,
 117            show arctan 0 - h_lower 0 = 0 by simp [h_lower, arctan_zero]]
 118
 119/-! ## §2. Constructive arctan(1/3) interval -/
 120
 121/-- Interval containing arctan(1/3) ≈ 0.32175 -/
 122def arctanOneThirdInterval : Interval where
 123  lo := 321 / 1000  -- 0.321
 124  hi := 322 / 1000  -- 0.322
 125  valid := by norm_num
 126
 127/-- **PROVED**: arctan(1/3) is in arctanOneThirdInterval.
 128    Lower bound from `h_lower(1/3) = 24628/76545 ≥ 321/1000`.
 129    Upper bound from `g_upper(1/3) = 391/1215 ≤ 322/1000`. -/
 130theorem arctan_one_third_in_interval :
 131    arctanOneThirdInterval.contains (arctan (1 / 3)) := by
 132  simp only [contains, arctanOneThirdInterval]
 133  constructor
 134  · -- lower: ↑(321/1000 : ℚ) ≤ arctan(1/3)
 135    have hlo := lower_poly_le_arctan (1 / 3) (by norm_num)
 136    have hval : h_lower (1 / 3) = 24628 / 76545 := by unfold h_lower; norm_num
 137    rw [hval] at hlo
 138    have hq : ((321 / 1000 : ℚ) : ℝ) ≤ (24628 : ℝ) / 76545 := by push_cast; norm_num
 139    linarith
 140  · -- upper: arctan(1/3) ≤ ↑(322/1000 : ℚ)
 141    have hhi := arctan_le_upper_poly (1 / 3) (by norm_num)
 142    have hval : g_upper (1 / 3) = 391 / 1215 := by unfold g_upper; norm_num
 143    rw [hval] at hhi
 144    have hq : (391 : ℝ) / 1215 ≤ ((322 / 1000 : ℚ) : ℝ) := by push_cast; norm_num
 145    linarith
 146
 147/-! ## §3. Constructive arctan(2) interval via addition formula -/
 148
 149/-- **PROVED**: `arctan(2) = π/4 + arctan(1/3)`.
 150    From the addition formula: `arctan(1) + arctan(1/3) = arctan(2)`. -/
 151theorem arctan_two_eq : arctan 2 = Real.pi / 4 + arctan (1 / 3) := by
 152  have h1 : arctan 1 = Real.pi / 4 := by
 153    rw [← tan_pi_div_four]
 154    exact arctan_tan (by linarith [pi_pos]) (by linarith [pi_pos])
 155  have hadd : arctan 1 + arctan (1 / 3) = arctan 2 := by
 156    rw [arctan_add (by norm_num : (1 : ℝ) * (1 / 3) < 1)]
 157    congr 1; norm_num
 158  linarith
 159
 160/-- Interval containing arctan 2 = π/4 + arctan(1/3) -/
 161def arctanTwoInterval : Interval :=
 162  let pi4 := smulPos (1 / 4) piInterval (by norm_num)
 163  add pi4 arctanOneThirdInterval
 164
 165/-- **PROVED**: arctan(2) is in arctanTwoInterval. -/
 166theorem arctan_two_in_interval :
 167    arctanTwoInterval.contains (arctan 2) := by
 168  rw [arctan_two_eq]
 169  unfold arctanTwoInterval
 170  apply add_contains_add
 171  · -- π/4 is in (1/4) · piInterval
 172    have hpi := pi_in_piInterval
 173    have := smulPos_contains_smul (q := 1 / 4) (by norm_num : (0 : ℚ) < 1 / 4) hpi
 174    simp only [Rat.cast_div, Rat.cast_one, Rat.cast_ofNat] at this
 175    convert this using 1
 176    ring
 177  · exact arctan_one_third_in_interval
 178
 179end IndisputableMonolith.Numerics
 180