import Omega.Affine
import Omega.Floor
import Omega.Ledger

/-! External scalar interface only. This file does not implement, modify, or
claim an extension of frozen SRSD v1.0. -/
namespace Omega
noncomputable section

def scalarEncode (h y : ℝ) : ℤ × ℝ := (Int.floor (y/h), y-h*Int.floor (y/h))
def scalarDecode (h : ℝ) (p : ℤ × ℝ) : ℝ := h*p.1+p.2

theorem scalar_roundtrip (h y : ℝ) : scalarDecode h (scalarEncode h y)=y := by
  simp [scalarEncode,scalarDecode]

theorem scalar_bounds (h y : ℝ) (hh : 0 < h) :
    0 ≤ (scalarEncode h y).2 ∧ (scalarEncode h y).2 < h := by
  have hl := Int.floor_le (y/h)
  have hu := Int.lt_floor_add_one (y/h)
  have hl' := (le_div_iff₀ hh).mp hl
  have hu' := (div_lt_iff₀ hh).mp hu
  dsimp [scalarEncode]
  constructor <;> nlinarith

theorem scalar_canonical (h η : ℝ) (e : ℤ) (hh : 0 < h)
    (hη : 0 ≤ η ∧ η < h) : scalarEncode h (h*e+η)=(e,η) := by
  have he : Int.floor ((h*(e:ℝ)+η)/h)=e := by
    apply Int.floor_eq_iff.mpr
    constructor
    · apply (le_div_iff₀ hh).mpr; nlinarith
    · apply (div_lt_iff₀ hh).mpr; nlinarith
  simp [scalarEncode,he]

theorem rational_residual (d : ℕ) : q*d-(F d:ℝ)=(rho d:ℝ)/S := by
  have hi : (S:ℝ)*(F d:ℝ)+(rho d:ℝ)=(Q:ℝ)*d := by
    exact_mod_cast division_identity d
  unfold q
  field_simp [show (S:ℝ)≠0 by norm_num [S]]
  nlinarith

theorem real_residual (d : ℕ) (hd : d < domainSize) :
    lam*d-(Int.floor (lam*d):ℝ)=(rho d:ℝ)/S+(lam-q)*d := by
  rw [certified_floor_equality d hd]
  push_cast
  rw [← rational_residual]
  ring

theorem positive_residual_difference (d : ℕ) (hd : 0 < d) :
    0 < (lam-q)*(d:ℝ) := mul_pos (sub_pos.mpr q_lt_lam) (by exact_mod_cast hd)

theorem grid_finite_interface (h : ℝ) (hh : 0 < h) (d : ℕ) (hd : d < domainSize) :
    quantizedStep h (h*d)=h*(F d:ℝ) := by
  have he := grid_interface h hh (d:ℤ)
  push_cast at he
  rw [he,certified_floor_equality d hd]
  norm_cast

theorem grid_trajectory (h : ℝ) (hh : 0 < h) (n d : ℕ) (hd : d < domainSize) :
    (quantizedStep h)^[n] (h*d)=h*(F^[n] d:ℝ) := by
  induction n generalizing d with
  | zero => rfl
  | succ n ih =>
    rw [Function.iterate_succ_apply,Function.iterate_succ_apply,
      grid_finite_interface h hh d hd]
    exact ih _ (forward_invariant d hd)

/-- A witness that projecting each step differs from rounding once at the end. -/
theorem repeated_rounding_counterexample :
    (realStep lam)^[2] 26=2 ∧ Int.floor (lam^2*26)=3 := by
  constructor
  · rw [trajectory_equality lam ⟨lam_in_interval.1.le,lam_in_interval.2⟩
      2 26 (by norm_num [domainSize])]
    norm_num [Function.iterate_succ_apply,F,Q,S]
  · apply Int.floor_eq_iff.mpr
    norm_num
    have hl : (0.3399:ℝ)<lam := by nlinarith [lam_bounds.1]
    have hu : lam<(0.34:ℝ) := by nlinarith [lam_bounds.2]
    constructor <;> nlinarith [lam_pos]

/-- The canonical Euclidean arithmetic is exactly the general ledger interface. -/
theorem ledger_interface (d : ℕ) :
    Ledger.step Q S d=(F d:ℤ) ∧ Ledger.remainder Q S d=(rho d:ℤ) := by
  simp [Ledger.step,Ledger.remainder,F,rho,Int.natCast_div,Int.natCast_mod]

theorem external_factorization (d : ℕ) :
    scalarEncode S ((Q:ℝ)*d)=((F d:ℤ),(rho d:ℝ)) := by
  have he : (Q:ℝ)*d=(S:ℝ)*(F d:ℤ)+(rho d:ℝ) := by
    push_cast
    exact_mod_cast (division_identity d).symm
  rw [he]
  apply scalar_canonical
  · norm_num [S]
  · constructor
    · positivity
    · exact_mod_cast remainder_bound d

theorem predecessor_inequalities (A B d e : ℤ) (hB : 0 < B) :
    Ledger.step A B d=e ↔ B*e ≤ A*d ∧ A*d < B*(e+1) := by
  rw [Ledger.step,Int.ediv_eq_iff_of_pos hB]
  constructor <;> intro h <;> constructor <;> nlinarith [h.1,h.2]

variable {V : Type*} [AddCommGroup V] [Module ℝ V]

/-- Reconstruction of points explicitly consumes the retained grid scale,
center and frame. The ledger alone reconstructs the integer state. -/
def pointFromState (h : ℝ) (α : V) (frame : ℝ →ₗ[ℝ] V) (d : ℤ) : V :=
  α+frame (h*d)

theorem point_reconstruction (h : ℝ) (α : V) (frame : ℝ →ₗ[ℝ] V)
    (n : ℕ) (d : ℤ) :
    pointFromState h α frame
      (Ledger.recover Q S (Ledger.encode Q S n d).1 (Ledger.encode Q S n d).2)=
      pointFromState h α frame d := by
  rw [Ledger.roundtrip (Q:ℤ) (S:ℤ) (by norm_num [Q])]

end
end Omega