import Mathlib.Data.Real.Pi.Bounds
import Mathlib.Tactic.Linarith
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Positivity
import Mathlib.Tactic.Ring

/-! Exact real constants. All decimal bounds below denote rational numbers.
No floating point, external evaluator, or project axioms are used. -/
namespace Omega
noncomputable section

def phi : ℝ := (1 + Real.sqrt 5) / 2
def epsilon : ℝ := Real.pi / 648000
def delta : ℝ := epsilon / 2
def s : ℝ := Real.sin delta
def c : ℝ := Real.cos delta
def lam : ℝ := phi / (phi + Real.pi)
def mu : ℝ := Real.pi / (phi + Real.pi)
def Q : ℕ := 339949771344778
def S : ℕ := 1000000000000000
def domainSize : ℕ := 32768
def q : ℝ := (Q : ℝ) / S

theorem phi_pos : 0 < phi := by unfold phi; positivity

theorem phi_sq : phi ^ 2 = phi + 1 := by
  have h := Real.sq_sqrt (show (0 : ℝ) ≤ 5 by norm_num)
  unfold phi
  nlinarith

theorem phi_bounds : 1.61803398874989484820 < phi ∧
    phi < 1.61803398874989484821 := by
  have h := phi_sq
  have hp := phi_pos
  constructor <;> nlinarith

theorem denominator_pos : 0 < phi + Real.pi := add_pos phi_pos Real.pi_pos

theorem lam_pos : 0 < lam := div_pos phi_pos denominator_pos

theorem lam_lt_one : lam < 1 := by
  unfold lam
  rw [div_lt_one denominator_pos]
  linarith [Real.pi_pos]

theorem shares_sum : lam + mu = 1 := by
  unfold lam mu
  rw [← add_div, div_self (ne_of_gt denominator_pos)]

theorem mu_eq_one_sub : mu = 1 - lam := by linarith [shares_sum]

theorem lam_bounds : (339949771344778 : ℝ) / 10^15 < lam ∧
    lam < (339949771344779 : ℝ) / 10^15 := by
  obtain ⟨hl, hu⟩ := phi_bounds
  have pl := Real.pi_gt_d20
  have pu := Real.pi_lt_d20
  unfold lam
  constructor
  · rw [lt_div_iff₀ denominator_pos]
    norm_num at *
    nlinarith
  · rw [div_lt_iff₀ denominator_pos]
    norm_num at *
    nlinarith

theorem q_lt_lam : q < lam := by
  have h := lam_bounds.1
  norm_num [q, Q, S] at *
  exact h

theorem lam_lt_q_plus : lam < q + 1 / S := by
  have h := lam_bounds.2
  norm_num [q, Q, S] at *
  exact h

theorem delta_pos : 0 < delta := by unfold delta epsilon; positivity

theorem delta_lt_half_pi : delta < Real.pi / 2 := by
  unfold delta epsilon
  linarith [Real.pi_pos]

theorem s_pos : 0 < s := by
  apply Real.sin_pos_of_pos_of_lt_pi delta_pos
  linarith [delta_lt_half_pi, Real.pi_pos]

theorem c_pos : 0 < c := by
  apply Real.cos_pos_of_mem_Ioo
  constructor
  · linarith [delta_pos, Real.pi_pos]
  · exact delta_lt_half_pi

theorem sc_identity : s ^ 2 + c ^ 2 = 1 := Real.sin_sq_add_cos_sq delta

theorem one_arcsecond : epsilon * (180 / Real.pi) * 3600 = 1 := by
  unfold epsilon
  field_simp
  <;> ring

end
end Omega