import Omega.Constants
import Mathlib.Algebra.Order.Floor.Ring
import Mathlib.Algebra.Order.Floor.Semiring

namespace Omega

def F (d : ℕ) : ℕ := Q*d/S
def rho (d : ℕ) : ℕ := Q*d%S

theorem S_pos : 0 < S := by norm_num [S]
theorem Q_pos : 0 < Q := by norm_num [Q]
theorem Q_lt_S : Q < S := by norm_num [Q,S]

theorem division_identity (d : ℕ) : S*F d + rho d = Q*d := by
  exact Nat.div_add_mod (Q*d) S

theorem remainder_bound (d : ℕ) : rho d < S := Nat.mod_lt _ S_pos

theorem step_zero : F 0 = 0 := by simp [F]

theorem step_lt (d : ℕ) (hd : 0 < d) : F d < d := by
  apply (Nat.div_lt_iff_lt_mul S_pos).2
  simpa [Nat.mul_comm] using Nat.mul_lt_mul_of_pos_right Q_lt_S hd

theorem step_le (d : ℕ) : F d ≤ d := by
  by_cases hd : d = 0
  · subst d; simp [step_zero]
  · exact Nat.le_of_lt (step_lt d (Nat.pos_of_ne_zero hd))

theorem forward_invariant (d : ℕ) (hd : d < domainSize) : F d < domainSize :=
  lt_of_le_of_lt (step_le d) hd

theorem step_monotone : Monotone F := by
  intro d e h
  exact Nat.div_le_div_right (Nat.mul_le_mul_left Q h)

noncomputable section

def lower : ℝ := 4196/12343
def upper : ℝ := 7039/20706

theorem q_in_interval : lower < q ∧ q < upper := by norm_num [lower,upper,q,Q,S]
theorem lam_in_interval : lower < lam ∧ lam < upper := by
  have := lam_bounds
  norm_num [lower,upper] at *
  constructor <;> linarith

theorem floor_q (d : ℕ) : Int.floor (q*d) = (F d : ℤ) := by
  apply Int.floor_eq_iff.mpr
  have he : (S:ℝ)*(F d:ℝ)+(rho d:ℝ) = (Q:ℝ)*d := by
    exact_mod_cast division_identity d
  have hb : (rho d:ℝ) < S := by exact_mod_cast remainder_bound d
  have hz : (0:ℝ) ≤ rho d := by positivity
  have hs : (0:ℝ) < S := by exact_mod_cast S_pos
  push_cast
  unfold q
  rw [div_mul_eq_mul_div]
  constructor
  · apply (le_div_iff₀ hs).2
    nlinarith
  · apply (div_lt_iff₀ hs).2
    nlinarith

/-- Farey determinant argument. No enumeration or external numerical certificate. -/
theorem no_interior_boundary (d : ℕ) (hd : d < domainSize) (k : ℤ)
    (hl : lower*d < k) (hu : (k:ℝ) < upper*d) : False := by
  have hxR : (0:ℝ) < 12343*(k:ℝ)-4196*(d:ℝ) := by
    unfold lower at hl
    linarith
  have hyR : (0:ℝ) < 7039*(d:ℝ)-20706*(k:ℝ) := by
    unfold upper at hu
    linarith
  have hx : (0:ℤ) < 12343*k-4196*d := by exact_mod_cast hxR
  have hy : (0:ℤ) < 7039*(d:ℤ)-20706*k := by exact_mod_cast hyR
  have hx1 : (1:ℤ) ≤ 12343*k-4196*d := by omega
  have hy1 : (1:ℤ) ≤ 7039*(d:ℤ)-20706*k := by omega
  have hi : (d:ℤ) = 20706*(12343*k-4196*d)+12343*(7039*d-20706*k) := by ring
  have hdI : (d:ℤ) < 32768 := by exact_mod_cast hd
  linarith

theorem floor_agreement_of_interval (r : ℝ) (hr : lower ≤ r ∧ r < upper)
    (d : ℕ) (hd : d < domainSize) : Int.floor (r*d) = (F d:ℤ) := by
  classical
  by_cases hz : d=0
  · subst d; simp [F]
  have hp : (0:ℝ) < d := by exact_mod_cast Nat.pos_of_ne_zero hz
  have hf := Int.floor_eq_iff.mp (floor_q d)
  push_cast at hf
  apply Int.floor_eq_iff.mpr
  push_cast
  constructor
  · by_contra h
    have hbad : r*(d:ℝ) < (F d:ℝ) := lt_of_not_ge h
    apply no_interior_boundary d hd (F d)
    · push_cast
      nlinarith [hr.1]
    · push_cast
      nlinarith [q_in_interval.2, hf.1]
  · by_contra h
    have hbad : (F d:ℝ)+1 ≤ r*d := le_of_not_gt h
    apply no_interior_boundary d hd ((F d:ℤ)+1)
    · push_cast
      nlinarith [q_in_interval.1, hf.2]
    · push_cast
      nlinarith [hr.2]

theorem exact_parameter_interval (r : ℝ) :
    (∀ d : ℕ, d < domainSize → Int.floor (r*d) = (F d:ℤ)) ↔
      lower ≤ r ∧ r < upper := by
  constructor
  · intro h
    have h1 := Int.floor_eq_iff.mp (h 12343 (by norm_num [domainSize]))
    have h2 := Int.floor_eq_iff.mp (h 20706 (by norm_num [domainSize]))
    norm_num [F,Q,S] at h1 h2
    constructor
    · unfold lower; linarith [h1.1]
    · unfold upper; linarith [h2.2]
  · exact fun hr d hd => floor_agreement_of_interval r hr d hd

theorem certified_floor_equality (d : ℕ) (hd : d < domainSize) :
    Int.floor (lam*d) = (F d:ℤ) :=
  floor_agreement_of_interval lam ⟨lam_in_interval.1.le,lam_in_interval.2⟩ d hd

def realStep (r : ℝ) (d : ℕ) : ℕ := (Int.floor (r*d)).toNat

theorem trajectory_equality (r : ℝ) (hr : lower ≤ r ∧ r < upper)
    (n d : ℕ) (hd : d < domainSize) : (realStep r)^[n] d = F^[n] d := by
  induction n generalizing d with
  | zero => rfl
  | succ n ih =>
    rw [Function.iterate_succ_apply, Function.iterate_succ_apply]
    have he : realStep r d = F d := by
      simp [realStep, floor_agreement_of_interval r hr d hd]
    rw [he]
    exact ih _ (forward_invariant d hd)

def truncation (p : ℕ) : ℝ := (Int.floor ((10:ℝ)^p*lam):ℝ) / (10:ℝ)^p

theorem truncation_le (p : ℕ) : truncation p ≤ lam := by
  unfold truncation
  apply (div_le_iff₀ (by positivity : (0:ℝ)<10^p)).2
  simpa [mul_comm] using Int.floor_le ((10:ℝ)^p*lam)

theorem truncation_succ (p : ℕ) : truncation p ≤ truncation (p+1) := by
  have hp : (0:ℝ)<10^p := by positivity
  have hf := Int.floor_le ((10:ℝ)^p*lam)
  have hi : 10*Int.floor ((10:ℝ)^p*lam) ≤ Int.floor ((10:ℝ)^(p+1)*lam) := by
    apply Int.le_floor.mpr
    push_cast
    rw [pow_succ]
    nlinarith
  have hc : (10:ℝ)*(Int.floor ((10:ℝ)^p*lam):ℝ) ≤
      (Int.floor ((10:ℝ)^(p+1)*lam):ℝ) := by exact_mod_cast hi
  unfold truncation
  apply (div_le_div_iff₀ hp (by positivity)).2
  rw [pow_succ] at hc ⊢
  nlinarith [mul_le_mul_of_nonneg_right hc hp.le]

theorem truncation_monotone : Monotone truncation := monotone_nat_of_le_succ truncation_succ

theorem truncation_seven : truncation 7 = 3399497/10000000 := by
  have h : Int.floor ((10:ℝ)^7*lam) = 3399497 := by
    apply Int.floor_eq_iff.mpr
    norm_num
    constructor <;> nlinarith [lam_bounds.1,lam_bounds.2]
  simp [truncation,h]
  <;> norm_num

theorem truncation_eight : truncation 8 = 33994977/100000000 := by
  have h : Int.floor ((10:ℝ)^8*lam) = 33994977 := by
    apply Int.floor_eq_iff.mpr
    norm_num
    constructor <;> nlinarith [lam_bounds.1,lam_bounds.2]
  simp [truncation,h]
  <;> norm_num

/-- All decimal precisions, rather than only the finite sweep through 18. -/
theorem precision_threshold (p : ℕ) :
    (∀ d : ℕ, d < domainSize → Int.floor (truncation p*d) = (F d:ℤ)) ↔ 8 ≤ p := by
  rw [exact_parameter_interval]
  constructor
  · intro h
    by_contra hp
    have hm := truncation_monotone (show p ≤ 7 by omega)
    rw [truncation_seven] at hm
    unfold lower at h
    linarith [h.1]
  · intro hp
    have hm := truncation_monotone hp
    rw [truncation_eight] at hm
    constructor
    · unfold lower
      linarith
    · exact lt_of_le_of_lt (truncation_le p) lam_in_interval.2

end
end Omega