SilverSight/formal/CoreFormalism/FixedPoint.lean

import Lean.Data.Json
import Mathlib.Data.UInt
import Mathlib.Tactic
import Mathlib.Data.Int.Basic
import Mathlib.Data.Nat.Basic

set_option maxRecDepth 20000
set_option linter.unusedSimpArgs false

namespace SilverSight.FixedPoint

open Lean

/-!
A proof-friendly fixed-point core.

Design rule:
* The semantic value is a bounded signed raw integer.
* Saturation is performed by `ofRawInt`.
* UInt bit-patterns are boundary/hardware artifacts, not the proof model.

This removes the old proof debt caused by proving signed arithmetic facts directly
against modular UInt32/UInt64 overflow behavior.
-/

-- ═══════════════════════════════════════════════════════════════════════════
-- Q0.16 signed normalized fraction
-- ═══════════════════════════════════════════════════════════════════════════

def q0_16MinRaw : Int := -32768
def q0_16MaxRaw : Int := 32767
def q0_16Scale  : Int := 32767

/--
Q0.16 pure fraction representation.
The canonical proof model stores the signed raw integer in [-32768, 32767].
Use boundary conversion functions when a UInt16 bit pattern is required.
-/
abbrev Q0_16 := { x : Int // q0_16MinRaw ≤ x ∧ x ≤ q0_16MaxRaw }

instance : Repr Q0_16 where
  reprPrec q _ := repr q.val

instance : BEq Q0_16 where
  beq a b := a.val == b.val

instance : Inhabited Q0_16 where
  default := ⟨0, by constructor <;> norm_num [q0_16MinRaw, q0_16MaxRaw]⟩

instance : ToJson Q0_16 where
  toJson q := Json.mkObj [("val", toJson q.val)]

namespace Q0_16

@[ext]
theorem ext {a b : Q0_16} (h : a.val = b.val) : a = b := Subtype.ext h

@[inline]
def toInt (q : Q0_16) : Int := q.val

@[inline]
def ofRawInt (raw : Int) : Q0_16 :=
  if hhi : raw > q0_16MaxRaw then
    ⟨q0_16MaxRaw, by constructor <;> norm_num [q0_16MinRaw, q0_16MaxRaw]⟩
  else if hlo : raw < q0_16MinRaw then
    ⟨q0_16MinRaw, by constructor <;> norm_num [q0_16MinRaw, q0_16MaxRaw]⟩
  else
    ⟨raw, by
      constructor
      · dsimp [q0_16MinRaw, q0_16MaxRaw] at *; omega
      · dsimp [q0_16MinRaw, q0_16MaxRaw] at *; omega⟩

instance : FromJson Q0_16 where
  fromJson? j := do
    let raw : Int ← fromJson? (← j.getObjVal? "val")
    pure (ofRawInt raw)

def zero : Q0_16 := ofRawInt 0
def one  : Q0_16 := ofRawInt q0_16MaxRaw
def half : Q0_16 := ofRawInt 16383

def neg (x : Q0_16) : Q0_16 := ofRawInt (-x.toInt)
def add (a b : Q0_16) : Q0_16 := ofRawInt (a.toInt + b.toInt)
def sub (a b : Q0_16) : Q0_16 := ofRawInt (a.toInt - b.toInt)
def mul (a b : Q0_16) : Q0_16 := ofRawInt ((a.toInt * b.toInt) / q0_16Scale)
def div (a b : Q0_16) : Q0_16 :=
  if b.toInt = 0 then zero else ofRawInt ((a.toInt * q0_16Scale) / b.toInt)
def abs (x : Q0_16) : Q0_16 := if x.toInt < 0 then neg x else x

instance : Add Q0_16 where add := add
instance : Sub Q0_16 where sub := sub
instance : Mul Q0_16 where mul := mul
instance : Div Q0_16 where div := div
instance : Neg Q0_16 where neg := neg

def lt (a b : Q0_16) : Bool := a.toInt < b.toInt
def le (a b : Q0_16) : Bool := a.toInt ≤ b.toInt
def gt (a b : Q0_16) : Bool := b.toInt < a.toInt
def ge (a b : Q0_16) : Bool := b.toInt ≤ a.toInt

def toFloat (q : Q0_16) : Float :=
  Float.ofInt q.toInt / 32767.0

def ofFloat (f : Float) : Q0_16 :=
  if f.isNaN then zero
  else if f ≥ 1.0 then one
  else if f ≤ -1.0 then neg one
  else if f < 0.0 then
    ofRawInt (-(Int.ofNat ((-f * 32767.0).round.toUInt16.toNat)))
  else
    ofRawInt (Int.ofNat ((f * 32767.0).round.toUInt16.toNat))

def log2 (q : Q0_16) : Q0_16 :=
  if q.toInt = 0 then zero
  else
    let f := toFloat q
    if f ≤ 0.0 then zero else ofFloat (Float.log2 f)

def min (a b : Q0_16) : Q0_16 :=
  if a.toInt ≤ b.toInt then a else b

end Q0_16

-- ═══════════════════════════════════════════════════════════════════════════
-- Q16.16 signed fixed-point
-- ═══════════════════════════════════════════════════════════════════════════

def q16MinRaw : Int := -2147483648
def q16MaxRaw : Int := 2147483647
def q16Scale  : Int := 65536

/-- Saturating clamp of a raw integer into [q16MinRaw, q16MaxRaw].
    This is the pure-Int kernel of `Q16_16.ofRawInt`; all monotonicity
    reasoning is proved once here and reused by higher lemmas. -/
def q16Clamp (i : Int) : Int :=
  if i > q16MaxRaw then q16MaxRaw
  else if i < q16MinRaw then q16MinRaw
  else i

/-- `q16Clamp` is monotone: a ≤ b → q16Clamp a ≤ q16Clamp b. -/
theorem q16Clamp_monotone (a b : Int) (h : a ≤ b) : q16Clamp a ≤ q16Clamp b := by
  unfold q16Clamp
  by_cases ha_hi : a > q16MaxRaw
  · by_cases hb_hi : b > q16MaxRaw
    · simp [ha_hi, hb_hi]
    · by_cases hb_lo : b < q16MinRaw
      · simp [ha_hi, hb_hi, hb_lo]; dsimp [q16MinRaw, q16MaxRaw] at *; omega
      · simp [ha_hi, hb_hi, hb_lo]; dsimp [q16MaxRaw] at *; omega
  · by_cases ha_lo : a < q16MinRaw
    · by_cases hb_hi : b > q16MaxRaw
      · simp [ha_hi, ha_lo, hb_hi]; dsimp [q16MinRaw, q16MaxRaw] at *; omega
      · by_cases hb_lo : b < q16MinRaw
        · simp [ha_hi, ha_lo, hb_hi, hb_lo]
        · simp [ha_hi, ha_lo, hb_hi, hb_lo]; dsimp [q16MinRaw] at *; omega
    · by_cases hb_hi : b > q16MaxRaw
      · simp [ha_hi, ha_lo, hb_hi]; dsimp [q16MaxRaw] at *; omega
      · by_cases hb_lo : b < q16MinRaw
        · simp [ha_hi, ha_lo, hb_hi, hb_lo]; dsimp [q16MinRaw] at *; omega
        · simp [ha_hi, ha_lo, hb_hi, hb_lo]; exact h

/-- `q16Clamp` is idempotent on in-range values. -/
theorem q16Clamp_id_of_inRange (i : Int) (hlo : q16MinRaw ≤ i) (hhi : i ≤ q16MaxRaw) :
    q16Clamp i = i := by
  unfold q16Clamp
  simp [show ¬ i > q16MaxRaw from by omega, show ¬ i < q16MinRaw from by omega]

lemma q16Clamp_lower (x : Int) : q16MinRaw ≤ q16Clamp x := by
  unfold q16Clamp q16MinRaw q16MaxRaw; split_ifs <;> omega

lemma q16Clamp_upper (x : Int) : q16Clamp x ≤ q16MaxRaw := by
  unfold q16Clamp q16MinRaw q16MaxRaw; split_ifs <;> omega

lemma q16Clamp_nonneg_of_nonneg {x : Int} (hx : 0 ≤ x) : 0 ≤ q16Clamp x := by
  unfold q16Clamp q16MinRaw q16MaxRaw
  split_ifs <;> omega

lemma q16Clamp_idem (x : Int) : q16Clamp (q16Clamp x) = q16Clamp x := by
  have h_upper := q16Clamp_upper x
  have h_lower := q16Clamp_lower x
  unfold q16Clamp q16MinRaw q16MaxRaw
  split_ifs <;> omega

/-- `q16Clamp` is 1-Lipschitz (non-expansive): clamped distance ≤ raw distance. -/
lemma q16Clamp_lipschitz (A B : Int) : |q16Clamp A - q16Clamp B| ≤ |A - B| := by
  by_cases h : A ≤ B
  · have hAB : A - B ≤ 0 := by omega
    have h_clamp : q16Clamp A - q16Clamp B ≤ 0 := by
      have h_mono : q16Clamp A ≤ q16Clamp B := q16Clamp_monotone A B h
      omega
    rw [abs_of_nonpos hAB, abs_of_nonpos h_clamp]
    unfold q16Clamp
    split <;> split <;> omega
  · have hBA : B ≤ A := by omega
    have hAB_pos : 0 ≤ A - B := by omega
    have h_clamp_pos : 0 ≤ q16Clamp A - q16Clamp B := by
      have h_mono : q16Clamp A ≥ q16Clamp B := q16Clamp_monotone B A hBA
      omega
    rw [abs_of_nonneg hAB_pos, abs_of_nonneg h_clamp_pos]
    unfold q16Clamp
    split <;> split <;> omega

/--
Q16.16 fixed-point representation.
The canonical proof model stores the signed raw integer in
[-2147483648, 2147483647].

Hardware/serialization UInt32 bit patterns should enter through `ofBits` and
leave through `toBits`. All semantic proofs use `toInt`.
-/
abbrev Q16_16 := { x : Int // q16MinRaw ≤ x ∧ x ≤ q16MaxRaw }

instance : Repr Q16_16 where
  reprPrec q _ := repr q.val

instance : BEq Q16_16 where
  beq a b := a.val == b.val

instance : Inhabited Q16_16 where
  default := ⟨0, by constructor <;> norm_num [q16MinRaw, q16MaxRaw]⟩

instance : ToJson Q16_16 where
  toJson q := Json.mkObj [("val", toJson q.val)]

namespace Q16_16

@[ext]
theorem ext {a b : Q16_16} (h : a.val = b.val) : a = b := Subtype.ext h

@[inline]
def toInt (q : Q16_16) : Int := q.val

@[inline]
def ofRawInt (raw : Int) : Q16_16 :=
  if hhi : raw > q16MaxRaw then
    ⟨q16MaxRaw, by constructor <;> norm_num [q16MinRaw, q16MaxRaw]⟩
  else if hlo : raw < q16MinRaw then
    ⟨q16MinRaw, by constructor <;> norm_num [q16MinRaw, q16MaxRaw]⟩
  else
    ⟨raw, by
      constructor
      · dsimp [q16MinRaw, q16MaxRaw] at *; omega
      · dsimp [q16MinRaw, q16MaxRaw] at *; omega⟩

instance : FromJson Q16_16 where
  fromJson? j := do
    let raw : Int ← fromJson? (← j.getObjVal? "val")
    pure (ofRawInt raw)

/-- Decode a UInt32 two's-complement hardware bit pattern into the signed model. -/
@[inline]
def ofBits (u : UInt32) : Q16_16 :=
  let n := u.toNat
  if n ≥ 2147483648 then ofRawInt ((n : Int) - 4294967296)
  else ofRawInt (n : Int)

/-- Encode the signed model as a UInt32 two's-complement hardware bit pattern. -/
@[inline]
def toBits (q : Q16_16) : UInt32 := UInt32.ofInt q.toInt

def zero : Q16_16 := ⟨0, by constructor <;> norm_num [q16MinRaw, q16MaxRaw]⟩
def one  : Q16_16 := ⟨q16Scale, by constructor <;> norm_num [q16MinRaw, q16MaxRaw, q16Scale]⟩
def negOne : Q16_16 := ⟨-q16Scale, by constructor <;> norm_num [q16MinRaw, q16MaxRaw, q16Scale]⟩
def epsilon : Q16_16 := ⟨1, by constructor <;> norm_num [q16MinRaw, q16MaxRaw]⟩
def two : Q16_16 := ⟨2 * q16Scale, by constructor <;> norm_num [q16MinRaw, q16MaxRaw, q16Scale]⟩
def maxVal : Q16_16 := ⟨q16MaxRaw, by constructor <;> norm_num [q16MinRaw, q16MaxRaw]⟩
def minVal : Q16_16 := ⟨q16MinRaw, by constructor <;> norm_num [q16MinRaw, q16MaxRaw]⟩

/-- Saturating infinity/illegal sentinel. If you need the old 0xFFFFFFFF bit
sentinel, use `ofRawInt (-1)` or `ofBits 0xFFFFFFFF` explicitly. -/
def infinity : Q16_16 := maxVal

def scale : Nat := 65536

def ofNat (n : Nat) : Q16_16 := ofRawInt ((n : Int) * q16Scale)

def satFromNat (n : Nat) : Q16_16 := ofRawInt ((n : Int) * q16Scale)

def ofRatio (num : Nat) (den : Nat) : Q16_16 :=
  if den = 0 then zero
  else ofRawInt (Int.ofNat (num * scale / den))

instance : OfNat Q16_16 n where
  ofNat := ofNat n

@[inline]
def ofInt (n : Int) : Q16_16 := ofRawInt (n * q16Scale)

/-- Saturating addition. -/
@[inline]
def add (a b : Q16_16) : Q16_16 := ofRawInt (a.toInt + b.toInt)

/-- Saturating subtraction. -/
@[inline]
def sub (a b : Q16_16) : Q16_16 := ofRawInt (a.toInt - b.toInt)

/-- Saturating Q16.16 multiplication: raw result is `(a*b)/65536`. -/
@[inline]
def mul (a b : Q16_16) : Q16_16 := ofRawInt ((a.toInt * b.toInt) / q16Scale)

/-- Saturating Q16.16 division: raw result is `(a*65536)/b`. -/
@[inline]
def div (a b : Q16_16) : Q16_16 :=
  if b.toInt = 0 then infinity else ofRawInt ((a.toInt * q16Scale) / b.toInt)

@[inline]
def neg (q : Q16_16) : Q16_16 := ofRawInt (-q.toInt)

@[inline]
def abs (q : Q16_16) : Q16_16 := if q.toInt < 0 then neg q else q

@[inline]
def ofFloat (f : Float) : Q16_16 :=
  if f.isNaN || f ≥ 32768.0 then infinity
  else if f ≤ -32768.0 then minVal
  else if f < 0.0 then
    ofRawInt (-(Int.ofNat ((-f * 65536.0).floor.toUInt32.toNat)))
  else
    ofRawInt (Int.ofNat ((f * 65536.0).floor.toUInt32.toNat))

/-- Integer square root via Newton's method. Returns floor(√n).
    Terminates in at most 64 iterations (fuel-bounded).
    Division-by-zero is impossible: `x ≥ 1` is a loop invariant
    because the initial value `n/2+1 ≥ 1` and every refinement `x'`
    is the average of positive integers. -/
private def intSqrt (n : Int) : Int :=
  if n ≤ 0 then 0
  else
    let rec loop (x : Int) (fuel : Nat) : Int :=
      match fuel with
      | 0     => x
      | f + 1 =>
        let x' := (x + n / x) / 2
        if x' ≥ x then x else loop x' f
    loop (n / 2 + 1) 64

/-- Q16.16 square root via integer Newton's method.
    Computes floor(√(q.raw × 65536)) which is the Q16.16
    representation of √(q.raw/65536). -/
@[inline]
def sqrt (q : Q16_16) : Q16_16 :=
  if q.toInt ≤ 0 then zero
  else ofRawInt (intSqrt (q.toInt * q16Scale))

@[inline]
def toFloat (q : Q16_16) : Float :=
  Float.ofInt q.toInt / 65536.0

/-- Natural logarithm approximation around 1.0. -/
def ln (q : Q16_16) : Q16_16 :=
  let x := q.toInt
  if x ≤ 0 then zero
  else
    let y := x - q16Scale
    let y2 := (y * y) / q16Scale
    let y3 := (y * y2) / q16Scale
    ofRawInt (y - y2 / 2 + y3 / 3)

def log2 (q : Q16_16) : Q16_16 :=
  let ln2 : Q16_16 := ofRawInt 45426
  div (ln q) ln2

def expNeg (x : Q16_16) : Q16_16 :=
  if x.toInt ≥ 0x00030000 then zero
  else if x.toInt ≥ 0x00020000 then ofRawInt 0x00004D29
  else if x.toInt ≥ 0x00010000 then ofRawInt 0x0000C5C0
  else ofRawInt 0x0001C5C0

/-- Q16.16 exponential via Taylor series eˣ ≈ Σ xⁿ/n! (n=0..6).
    Valid for |x| ≤ 1. For larger x, range-reduce via eˣ = (e^(x/2))²
    (repeated squaring). Achieves ~Q16.16 precision. -/
def exp (x : Q16_16) : Q16_16 :=
  if x.toInt ≤ -q16Scale then zero
  else if x.toInt ≥ 4 * q16Scale then maxVal
  else
    -- Range-reduce: x = k*ln(2) + r, |r| < ln(2), then eˣ = 2ᵏ * eʳ
    let ln2Raw := 45426  -- Q16.16 representation of ln(2) ≈ 0.693147
    let k := x.toInt / ln2Raw
    let r := x.toInt - k * ln2Raw
    -- Taylor for eʳ through r⁶/720 (|r| < 0.693, error < 2e-6)
    let r2 := (r * r) / q16Scale
    let r3 := (r * r2) / q16Scale
    let r4 := (r2 * r2) / q16Scale
    let r5 := (r2 * r3) / q16Scale
    let r6 := (r3 * r3) / q16Scale
    let taylor := q16Scale + r + r2 / 2 + r3 / 6 + r4 / 24 + r5 / 120 + r6 / 720
    -- Scale by 2ᵏ (left-shift by k, saturate)
    let shifted := if k ≥ 15 then q16MaxRaw else if k ≤ -15 then 0
                   else if k ≥ 0 then taylor <<< k.toNat
                   else taylor >>> (-k).toNat
    ofRawInt (if shifted > q16MaxRaw then q16MaxRaw else if shifted < q16MinRaw then q16MinRaw else shifted)

/-- Q16.16 sine via 7th-order Taylor: sin(x) ≈ x - x³/6 + x⁵/120 - x⁷/5040.
    Range-reduced to [0, π/2] using symmetry. -/
def sin (x : Q16_16) : Q16_16 :=
  -- Map [−π, π] → [0, π], then reflect
  let piRaw := 205887  -- Q16.16 π ≈ 3.14159
  let twoPiRaw := 411774
  -- Normalize to [0, 2π)
  let raw := x.toInt % twoPiRaw
  let raw' := if raw < 0 then raw + twoPiRaw else raw
  -- Quadrant: 0 = [0,π/2), 1 = [π/2,π), 2 = [π,3π/2), 3 = [3π/2,2π)
  let halfPi := piRaw / 2
  let q := if raw' < halfPi then 0
           else if raw' < piRaw then 1
           else if raw' < piRaw + halfPi then 2
           else 3
  -- Reduce to [0, π/2]
  let reduced :=
    match q with
    | 0 => raw'
    | 1 => piRaw - raw'
    | 2 => raw' - piRaw
    | 3 => twoPiRaw - raw'
    | _ => raw'  -- unreachable
  -- Taylor sin(t) ≈ t - t³/6 + t⁵/120 - t⁷/5040 (|t| ≤ π/2 ≈ 1.57)
  let t := reduced
  let t2 := (t * t) / q16Scale
  let t3 := (t * t2) / q16Scale
  let t5 := (t3 * t2) / q16Scale
  let t7 := (t5 * t2) / q16Scale
  let sinPos := t - t3 / 6 + t5 / 120 - t7 / 5040
  -- Sign by quadrant: quadrants 0,1 → positive; quadrants 2,3 → negative
  if q < 2 then ofRawInt (max q16MinRaw (min q16MaxRaw sinPos))
  else ofRawInt (max q16MinRaw (min q16MaxRaw (-sinPos)))

/-- Q16.16 power: base^e = exp(e * ln(base)). -/
def pow (base e : Q16_16) : Q16_16 :=
  if base.toInt ≤ 0 then
    if e.toInt = 0 then one else zero
  else if e.toInt = 0 then one
  else if e.toInt = q16Scale then base  -- e = 1
  else exp (mul (ln base) e)

/-- Natural logarithm — alias for `ln`. -/
def log (x : Q16_16) : Q16_16 := ln x

-- ═══════════════════════════════════════════════════════════════════════════
-- Inverse Trigonometric Functions (integer-only, no Float)
-- ═══════════════════════════════════════════════════════════════════════════

/-- Core arctangent for |x| ≤ 1 using minimax polynomial.
    atan(x) ≈ x * (0.9998660 + x²*(-0.3302995 + x²*(0.1801410 + x²*(-0.0851330))))
    Max error < 2e-5 on [0,1]. -/
private def atanCore (x : Q16_16) : Q16_16 :=
  let xRaw := x.toInt
  -- x² in Q16.16
  let x2Raw := (xRaw * xRaw) / q16Scale
  -- Horner's method with step-by-step scaling to avoid overflow:
  --   p(x) = c0 + x²*(c1 + x²*(c2 + x²*c3))
  --   c0=65527, c1=-21642, c2=11804, c3=-5579
  let p3 : Int := -5579
  let p2 : Int := 11804 + (x2Raw * p3) / q16Scale
  let p1 : Int := -21642 + (x2Raw * p2) / q16Scale
  let p0 : Int := 65527 + (x2Raw * p1) / q16Scale
  ofRawInt ((xRaw * p0) / q16Scale)

/-- Q16.16 arctangent via minimax polynomial with range reduction.
    For |x| ≤ 1: uses atanCore directly.
    For |x| > 1: uses identity atan(x) = π/2 - atan(1/x).
    For x < 0: uses identity atan(x) = -atan(-x). -/
def atan (x : Q16_16) : Q16_16 :=
  if x.toInt = 0 then zero
  else
    let ax := abs x
    let halfPiRaw := 102943  -- Q16.16 π/2 ≈ 1.5708
    -- Range-reduce: if |x| > 1, use atan(x) = π/2 - atan(1/x)
    let result :=
      if ax.toInt ≤ q16Scale then
        atanCore ax
      else
        -- 1/x in Q16.16: (q16Scale² / ax_raw)
        let invX := ofRawInt ((q16Scale * q16Scale) / ax.toInt)
        sub (ofRawInt halfPiRaw) (atanCore invX)
    -- Sign: atan(-x) = -atan(x)
    if x.toInt < 0 then neg result else result

/-- Q16.16 arcsine via identity: asin(x) = atan(x / sqrt(1 - x²)).
    Valid for |x| ≤ 1. For |x| > 1, clips to ±π/2. -/
def asin (x : Q16_16) : Q16_16 :=
  let xRaw := x.toInt
  if xRaw ≥ q16Scale then div (ofRawInt 205887) two  -- π/2
  else if xRaw ≤ -q16Scale then neg (div (ofRawInt 205887) two)  -- -π/2
  else if xRaw = 0 then zero
  else
    -- Compute x / sqrt(1 - x²) in Q16.16
    -- 1 - x² in Q16.16: q16Scale - (xRaw*xRaw)/q16Scale
    let oneMinusX2 := q16Scale - (xRaw * xRaw) / q16Scale
    if oneMinusX2 ≤ 0 then
      -- Numerical underflow: |x| ≈ 1, return ±π/2
      if xRaw > 0 then div (ofRawInt 205887) two else neg (div (ofRawInt 205887) two)
    else
      let sqrtVal := sqrt (ofRawInt oneMinusX2)
      if sqrtVal.toInt = 0 then
        if xRaw > 0 then div (ofRawInt 205887) two else neg (div (ofRawInt 205887) two)
      else
        atan (div x sqrtVal)

/-- Q16.16 arccosine via identity: acos(x) = π/2 - asin(x).
    Valid for |x| ≤ 1. -/
def acos (x : Q16_16) : Q16_16 :=
  sub (div (ofRawInt 205887) two) (asin x)

/-- Q16.16 two-argument arctangent.
    Computes the angle (radians) from the positive x-axis to the point (x, y).
    Handles all four quadrants correctly. -/
def atan2 (y x : Q16_16) : Q16_16 :=
  let piRaw : Q16_16 := ofRawInt 205887
  if x.toInt = 0 then
    if y.toInt = 0 then zero  -- undefined, return 0
    else if y.toInt > 0 then div piRaw two  -- π/2
    else neg (div piRaw two)  -- -π/2
  else if x.toInt > 0 then
    atan (div y x)
  else if y.toInt ≥ 0 then
    add (atan (div y x)) piRaw  -- quadrant II: atan(y/x) + π
  else
    sub (atan (div y x)) piRaw  -- quadrant III: atan(y/x) - π

instance : Add Q16_16 := ⟨add⟩
instance : Sub Q16_16 := ⟨sub⟩
instance : Mul Q16_16 := ⟨mul⟩
instance : Div Q16_16 := ⟨div⟩
instance : Neg Q16_16 := ⟨neg⟩

instance : LE Q16_16 where
  le a b := a.toInt ≤ b.toInt

instance : LT Q16_16 where
  lt a b := a.toInt < b.toInt

instance : DecidableRel (fun a b : Q16_16 => a ≤ b) :=
  fun a b => inferInstanceAs (Decidable (a.toInt ≤ b.toInt))

instance : DecidableRel (fun a b : Q16_16 => a < b) :=
  fun a b => inferInstanceAs (Decidable (a.toInt < b.toInt))

@[inline]
def ge (a b : Q16_16) : Bool := b.toInt ≤ a.toInt

@[inline]
def gt (a b : Q16_16) : Bool := b.toInt < a.toInt

def lt (a b : Q16_16) : Bool := a.toInt < b.toInt

def le (a b : Q16_16) : Bool := a.toInt ≤ b.toInt

def isNeg (q : Q16_16) : Bool := q.toInt < 0

def clip (x lo hi : Q16_16) : Q16_16 :=
  if x.toInt < lo.toInt then lo
  else if x.toInt > hi.toInt then hi
  else x

def sat01 (q : Q16_16) : Q16_16 :=
  if q.toInt < 0 then zero
  else if q.toInt > q16Scale then one
  else q

def max (a b : Q16_16) : Q16_16 :=
  if a.toInt ≥ b.toInt then a else b

def min (a b : Q16_16) : Q16_16 :=
  if a.toInt ≤ b.toInt then a else b

def recip (x : Q16_16) : Q16_16 :=
  let xInt := x.toInt
  if xInt = 0 then maxVal
  else
    let numer : Int := 0x100000000
    let denom := if xInt < 0 then -xInt else xInt
    let r := numer / denom
    let y := ofRawInt r
    if xInt < 0 then neg y else y

def ofRaw (n : Nat) : Q16_16 := ofRawInt (n : Int)

-- ═══════════════════════════════════════════════════════════════════════════
-- Algebraic lemmas
-- ═══════════════════════════════════════════════════════════════════════════

@[simp] theorem zero_toInt : toInt zero = 0 := rfl
@[simp] theorem one_toInt : toInt one = 65536 := rfl
@[simp] theorem epsilon_toInt : toInt epsilon = 1 := rfl

theorem epsilon_toInt_pos : toInt epsilon > 0 := by norm_num [epsilon_toInt]

@[simp]
theorem maxVal_toInt : toInt maxVal = q16MaxRaw := rfl

@[simp]
theorem minVal_toInt : toInt minVal = q16MinRaw := rfl

@[simp]
private theorem ofRawInt_zero : ofRawInt 0 = zero := by
  apply Subtype.ext
  simp [ofRawInt, zero, toInt, q16MinRaw, q16MaxRaw]

/-- Saturation lower-bound preservation. -/
theorem ofRawInt_toInt_ge (i c : Int)
    (hi : i ≥ c) (hcMin : q16MinRaw ≤ c) (hcMax : c ≤ q16MaxRaw) :
    (ofRawInt i).toInt ≥ c := by
  unfold ofRawInt toInt
  by_cases hhi : i > q16MaxRaw
  · simp [hhi]
    dsimp [q16MaxRaw] at *
    omega
  · by_cases hlo : i < q16MinRaw
    · dsimp [q16MinRaw, q16MaxRaw] at *
      omega
    · simp [hhi, hlo]
      exact hi

/-- Saturation preserves non-negativity for non-negative raw values. -/
theorem ofRawInt_toInt_nonneg (i : Int) (hi : i ≥ 0) :
    (ofRawInt i).toInt ≥ 0 := by
  exact ofRawInt_toInt_ge i 0 hi (by norm_num [q16MinRaw]) (by norm_num [q16MaxRaw])

/-- Bounded raw reconstruction. -/
theorem ofRawInt_toInt (a : Q16_16) : ofRawInt a.toInt = a := by
  apply Subtype.ext
  unfold ofRawInt toInt
  have hhi : ¬ a.val > q16MaxRaw := by exact not_lt.mpr a.property.2
  have hlo : ¬ a.val < q16MinRaw := by exact not_lt.mpr a.property.1
  simp [hhi, hlo]

theorem ofRawInt_toInt_eq_nonneg (i : Int) (h1 : i ≥ 0) (h2 : i ≤ q16MaxRaw) :
    (ofRawInt i).toInt = i := by
  unfold ofRawInt toInt
  have hhi : ¬ i > q16MaxRaw := by omega
  have hlo : ¬ i < q16MinRaw := by
    dsimp [q16MinRaw] at *
    omega
  simp [hhi, hlo]

/-- Saturating constructor lemma for non-negative raw values. -/
theorem ofRawInt_toInt_eq_general (i : Int) (h1 : i ≥ 0) :
    (ofRawInt i).toInt = if i > q16MaxRaw then q16MaxRaw else i := by
  by_cases h : i > q16MaxRaw
  · unfold ofRawInt toInt
    simp [h]
  · have hlo : ¬ i < q16MinRaw := by
      dsimp [q16MinRaw, q16MaxRaw] at *
      omega
    unfold ofRawInt toInt
    simp [h, hlo]

/-- `(ofRawInt i).toInt = q16Clamp i` — the bridge between the subtype
    constructor and the pure-Int clamp function. -/
theorem ofRawInt_toInt_eq_clamp (i : Int) : (ofRawInt i).toInt = q16Clamp i := by
  unfold ofRawInt toInt q16Clamp
  split_ifs <;> rfl

/-- `@[simp]` version rewriting `.val` directly (avoids `toInt` unfolding
    ordering issues in `simp` calls). -/
@[simp] theorem ofRawInt_val_eq_q16Clamp (i : Int) : (ofRawInt i).val = q16Clamp i :=
  ofRawInt_toInt_eq_clamp i

/-- `ofRawInt` is monotone: a ≤ b → (ofRawInt a).toInt ≤ (ofRawInt b).toInt.
    One-liner via q16Clamp_monotone. -/
theorem ofRawInt_monotone (a b : Int) (h : a ≤ b) :
    (ofRawInt a).toInt ≤ (ofRawInt b).toInt := by
  simp only [ofRawInt_toInt_eq_clamp]
  exact q16Clamp_monotone a b h

/-- Adding a nonnegative Q16_16 value cannot decrease the saturated result.
    This is the general form of the motif-scoring monotonicity used in
    Semantics.PIST.Motif §6.2: motifScore(match=true) ≥ motifScore(match=false). -/
theorem add_nonneg_monotone (a b : Q16_16) (hb : 0 ≤ b.toInt) :
    a.toInt ≤ (add a b).toInt := by
  unfold add
  have := ofRawInt_monotone a.toInt (a.toInt + b.toInt) (by omega)
  rwa [ofRawInt_toInt] at this

/-- zero * a = zero. -/
theorem zero_mul (a : Q16_16) : mul zero a = zero := by
  unfold mul
  rw [zero_toInt]
  simp

/-- a * zero = zero. -/
theorem mul_zero (a : Q16_16) : mul a zero = zero := by
  unfold mul
  rw [zero_toInt]
  simp

/-- a - a = zero. -/
theorem sub_self (a : Q16_16) : sub a a = zero := by
  unfold sub
  simp

/-- a + zero = a. -/
theorem add_zero (a : Q16_16) : add a zero = a := by
  unfold add
  rw [zero_toInt]
  simp
  exact ofRawInt_toInt a

/-- zero + a = a. -/
theorem zero_add (a : Q16_16) : add zero a = a := by
  unfold add
  rw [zero_toInt]
  simp
  exact ofRawInt_toInt a

/-- sqrt zero is zero. -/
theorem sqrt_zero : sqrt zero = zero := by
  unfold sqrt
  simp

/-- sqrt one is within one LSB of one. -/
theorem sqrt_one : (sqrt one).toInt - one.toInt ≤ 1 := by
  native_decide

#eval sqrt zero            -- should be zero
#eval sqrt one             -- should be ~1.0 (LSB-aligned)
#eval sqrt (Q16_16.ofNat 4)  -- should be ~2.0
#eval sqrt (Q16_16.ofNat 9)  -- should be ~3.0
#eval (sqrt (Q16_16.ofNat 4)).toFloat ≤ 2.1   -- Boolean guard witness

private theorem int_scale_mul_ediv_cancel (n : Int) : (q16Scale * n) / q16Scale = n := by
  rw [Int.mul_ediv_cancel_left]
  norm_num [q16Scale]

/-- one * a = a. -/
theorem one_mul (a : Q16_16) : mul one a = a := by
  unfold mul
  show ofRawInt (one.toInt * a.toInt / q16Scale) = a
  rw [show one.toInt = q16Scale from rfl]
  have h : (q16Scale * a.toInt) / q16Scale = a.toInt := int_scale_mul_ediv_cancel a.toInt
  rw [h]
  exact ofRawInt_toInt a

/-- a * one = a. -/
theorem mul_one (a : Q16_16) : mul a one = a := by
  unfold mul
  show ofRawInt (a.toInt * one.toInt / q16Scale) = a
  rw [show one.toInt = q16Scale from rfl]
  have h : (a.toInt * q16Scale) / q16Scale = a.toInt := by
    rw [Int.mul_comm]
    exact int_scale_mul_ediv_cancel a.toInt
  rw [h]
  exact ofRawInt_toInt a

/-- toInt = 0 iff the value is zero. -/
theorem toInt_eq_zero_iff {a : Q16_16} : a.toInt = 0 ↔ a = zero := by
  constructor
  · intro h
    apply Subtype.ext
    simpa [toInt, zero] using h
  · intro h
    rw [h]
    rfl

/-- zero / x = zero for any nonzero denominator. -/
theorem zero_div (x : Q16_16) (hx : x.val ≠ 0) : div zero x = zero := by
  unfold div
  have hx' : ¬ x.toInt = 0 := by
    simpa [toInt] using hx
  simp [hx', zero_toInt]

/-- Square is non-negative under signed saturating multiplication. -/
theorem mul_self_nonneg (a : Q16_16) : (mul a a).toInt ≥ 0 := by
  unfold mul
  have hprod : a.toInt * a.toInt ≥ 0 := by nlinarith
  have hdiv : (a.toInt * a.toInt) / q16Scale ≥ 0 := by
    apply Int.ediv_nonneg
    · exact hprod
    · norm_num [q16Scale]
  exact ofRawInt_toInt_nonneg ((a.toInt * a.toInt) / q16Scale) hdiv

/-- Product of two non-negative Q16.16 values is non-negative. -/
theorem mul_toInt_nonneg (a b : Q16_16) (ha : a.toInt ≥ 0) (hb : b.toInt ≥ 0) :
    (mul a b).toInt ≥ 0 := by
  unfold mul
  have hprod : a.toInt * b.toInt ≥ 0 := by nlinarith
  have hdiv : (a.toInt * b.toInt) / q16Scale ≥ 0 := by
    apply Int.ediv_nonneg
    · exact hprod
    · norm_num [q16Scale]
  exact ofRawInt_toInt_nonneg ((a.toInt * b.toInt) / q16Scale) hdiv

/-- Non-negative addition stays non-negative under saturation. -/
theorem ofRaw_toInt_nonneg (acc wcc : Q16_16)
    (hacc : acc.toInt ≥ 0) (hwcc : wcc.toInt ≥ 0) :
    (Q16_16.add acc wcc).toInt ≥ 0 := by
  unfold add
  have hsum : acc.toInt + wcc.toInt ≥ 0 := by omega
  exact ofRawInt_toInt_nonneg (acc.toInt + wcc.toInt) hsum

/-- Compatibility lemma: a non-negative raw value decoded through saturation is non-negative. -/
theorem mk_lt_half_nonneg (s : Int) (hs : s ≥ 0) (_h : s < 2147483648) :
    (ofRawInt s).toInt ≥ 0 := by
  exact ofRawInt_toInt_nonneg s hs

/-- Positive raw addition remains positive under saturation. -/
theorem add_pos_of_pos (a b : Q16_16) (ha : a.toInt > 0) (hb : b.toInt > 0) :
    (add a b).toInt > 0 := by
  unfold add
  have hsum : a.toInt + b.toInt ≥ 1 := by omega
  have hge : (ofRawInt (a.toInt + b.toInt)).toInt ≥ 1 :=
    ofRawInt_toInt_ge (a.toInt + b.toInt) 1 hsum
      (by norm_num [q16MinRaw]) (by norm_num [q16MaxRaw])
  omega

/-- (1 + omega).toInt ≥ 65536 when omega.toInt ≥ 0. -/
theorem add_one_omega_ge_one (omega : Q16_16) (h : omega.toInt ≥ 0) :
    (add one omega).toInt ≥ 65536 := by
  unfold add
  have hsum : one.toInt + omega.toInt ≥ 65536 := by
    rw [one_toInt]
    omega
  exact ofRawInt_toInt_ge (one.toInt + omega.toInt) 65536 hsum
    (by norm_num [q16MinRaw]) (by norm_num [q16MaxRaw])

/-- Non-negative Q16.16 values are bounded by maxVal. -/
theorem toInt_nonneg_le_maxVal (q : Q16_16) (_h : q.toInt ≥ 0) : q.toInt ≤ q16MaxRaw := by
  exact q.property.2

/-- Adding epsilon to a non-negative value yields a positive value. -/
theorem epsilon_add_pos {r : Q16_16} (hr : r.toInt ≥ 0) :
    (r + epsilon).toInt > 0 := by
  change (add r epsilon).toInt > 0
  unfold add
  have hsum : r.toInt + epsilon.toInt ≥ 1 := by
    rw [epsilon_toInt]
    omega
  have hge : (ofRawInt (r.toInt + epsilon.toInt)).toInt ≥ 1 :=
    ofRawInt_toInt_ge (r.toInt + epsilon.toInt) 1 hsum
      (by norm_num [q16MinRaw]) (by norm_num [q16MaxRaw])
  omega

/-- `abs (sub a b) = abs (sub b a)` — absolute value of a difference
    is symmetric. Holds for all Q16_16 values (proved by case analysis
    on `a.val - b.val` at the Int clamping boundary). -/
lemma q16Clamp_eq_q16MaxRaw_of_ge {x : Int} (h : x ≥ q16MaxRaw) : q16Clamp x = q16MaxRaw := by
  dsimp [q16Clamp]
  by_cases hx : x > q16MaxRaw
  · simp [hx]
  · have hx_eq : x = q16MaxRaw := le_antisymm (le_of_not_gt hx) h
    subst hx_eq; simp [q16Clamp, q16MaxRaw, q16MinRaw]

lemma q16Clamp_eq_q16MinRaw_of_le {x : Int} (h : x ≤ q16MinRaw) : q16Clamp x = q16MinRaw := by
  dsimp [q16Clamp]
  by_cases hx_hi : x > q16MaxRaw
  · unfold q16MinRaw q16MaxRaw at *; omega
  · by_cases hx_lo : x < q16MinRaw
    · simp [hx_hi, hx_lo]
    · have hx_eq : x = q16MinRaw := le_antisymm h (by
        by_contra hlt
        apply hx_lo
        exact lt_of_not_ge hlt)
      subst hx_eq; simp [q16Clamp, q16MaxRaw, q16MinRaw]

private lemma not_q16MaxRaw_lt_0 : ¬ q16MaxRaw < 0 := by unfold q16MaxRaw; omega

private lemma q16Clamp_2147483648_eq_q16MaxRaw : q16Clamp (2147483648 : Int) = q16MaxRaw := by
  unfold q16Clamp q16MaxRaw q16MinRaw; omega

private lemma val_if (c : Prop) [Decidable c] (x y : Q16_16) : (if c then x else y).val = (if c then x.val else y.val) := by
  split <;> rfl

theorem abs_sub_comm (a b : Q16_16) : abs (sub a b) = abs (sub b a) := by
  apply Q16_16.ext
  simp [abs, sub, neg, toInt, val_if, ofRawInt_val_eq_q16Clamp]
  have hswap : q16Clamp (b.val - a.val) = q16Clamp (-(a.val - b.val)) := by
    have : b.val - a.val = -(a.val - b.val) := by omega
    rw [this]
  rw [hswap]
  set d := a.val - b.val
  have hswap' : q16Clamp (-(a.val - b.val)) = q16Clamp (-d) := by rfl
  rw [hswap']
  by_cases hd_above : d > q16MaxRaw
  · have h_cd : q16Clamp d = q16MaxRaw := q16Clamp_eq_q16MaxRaw_of_ge (le_of_lt hd_above)
    have h_nd_low : -d ≤ q16MinRaw := by unfold q16MaxRaw q16MinRaw at *; omega
    have h_cnd : q16Clamp (-d) = q16MinRaw := q16Clamp_eq_q16MinRaw_of_le h_nd_low
    rw [h_cd, h_cnd]
    unfold q16Clamp q16MaxRaw q16MinRaw; norm_num
  · by_cases hd_below : d < q16MinRaw
    · have h_cd : q16Clamp d = q16MinRaw := q16Clamp_eq_q16MinRaw_of_le (by omega)
      have h_nd_high : -d ≥ q16MaxRaw := by unfold q16MaxRaw q16MinRaw at *; omega
      have h_cnd : q16Clamp (-d) = q16MaxRaw := q16Clamp_eq_q16MaxRaw_of_ge h_nd_high
      rw [h_cd, h_cnd]
      unfold q16Clamp q16MaxRaw q16MinRaw; norm_num
    · have h_lo : q16MinRaw ≤ d := by unfold q16MinRaw q16MaxRaw at *; omega
      have h_hi : d ≤ q16MaxRaw := by unfold q16MinRaw q16MaxRaw at *; omega
      have h_cd : q16Clamp d = d := q16Clamp_id_of_inRange d h_lo h_hi
      rw [h_cd]
      by_cases h_nd_above : -d > q16MaxRaw
      · have h_d_min : d = q16MinRaw := by unfold q16MinRaw q16MaxRaw at *; omega
        rw [h_d_min]
        unfold q16Clamp q16MaxRaw q16MinRaw; norm_num
      · by_cases h_nd_below : -d < q16MinRaw
        · have h_d_max : d = q16MaxRaw := by unfold q16MinRaw q16MaxRaw at *; omega
          rw [h_d_max]
          unfold q16Clamp q16MaxRaw q16MinRaw; norm_num
        · have h_nd_lo' : q16MinRaw ≤ -d := by omega
          have h_nd_hi' : -d ≤ q16MaxRaw := by omega
          have h_cnd : q16Clamp (-d) = -d := q16Clamp_id_of_inRange (-d) h_nd_lo' h_nd_hi'
          rw [h_cnd, show (-(-d : Int) = d) by omega]
          by_cases hd_neg : d < 0
          · omega
          · by_cases hd_zero : d = 0
            · rw [hd_zero]
              unfold q16Clamp q16MinRaw q16MaxRaw; norm_num
            · have hd_pos : 0 < d := by omega
              rw [h_cd]
              split_ifs <;> omega

/-- Subtraction is addition of the negation: a - b = a + (-b).

    NOTE: This theorem is NOT universally true for Q16_16. Counterexample:
    `a = b = q16MinRaw` gives LHS = 0, RHS = -1. The difference arises because
    `neg q16MinRaw` overflows to `q16MaxRaw`, altering the clamping path.
    SSMS does not use this theorem — the `bound` proof has been restructured
    to avoid it. -/
theorem sub_eq_add_neg (a b : Q16_16) (hb : b.toInt > q16MinRaw) : sub a b = add a (neg b) := by
  have h_neg_inRange_lo : q16MinRaw ≤ -b.toInt := by
    have h := b.property.2
    dsimp [toInt] at h ⊢
    dsimp [q16MaxRaw, q16MinRaw] at h ⊢
    omega
  have h_neg_inRange_hi : -b.toInt ≤ q16MaxRaw := by
    dsimp [toInt] at hb ⊢
    dsimp [q16MinRaw] at hb
    dsimp [q16MaxRaw]
    omega
  have h_neg_int : (neg b).toInt = -b.toInt := by
    rw [neg, ofRawInt_toInt_eq_clamp]
    apply q16Clamp_id_of_inRange _ h_neg_inRange_lo h_neg_inRange_hi
  rw [sub, add, h_neg_int]
  rfl

/-- Multiplication by a non-negative scalar is monotone:
    if a ≤ b and c ≥ 0, then a*c ≤ b*c.
    Used in SSMS.aciPreservedByMlgruStep for bound propagation. -/
theorem mul_mono_left (a b c : Q16_16) (h : a.toInt ≤ b.toInt) (hc : c.toInt ≥ 0) :
    (mul a c).toInt ≤ (mul b c).toInt := by
  unfold mul
  have hmul : a.toInt * c.toInt ≤ b.toInt * c.toInt := by
    apply Int.mul_le_mul_of_nonneg_right h hc
  have hpos : 0 < q16Scale := by unfold q16Scale; norm_num
  have hdiv : (a.toInt * c.toInt) / q16Scale ≤ (b.toInt * c.toInt) / q16Scale := by
    apply Int.ediv_le_ediv hpos hmul
  rw [ofRawInt_toInt_eq_clamp, ofRawInt_toInt_eq_clamp]
  exact q16Clamp_monotone _ _ hdiv

/-- Multiplication by a non-negative scalar is monotone on the right. -/
theorem mul_mono_right (a b c : Q16_16) (h : a.toInt ≤ b.toInt) (hc : c.toInt ≥ 0) :
    (mul c a).toInt ≤ (mul c b).toInt := by
  unfold mul
  have hmul : c.toInt * a.toInt ≤ c.toInt * b.toInt := by
    apply Int.mul_le_mul_of_nonneg_left h hc
  have hpos : 0 < q16Scale := by unfold q16Scale; norm_num
  have hdiv : (c.toInt * a.toInt) / q16Scale ≤ (c.toInt * b.toInt) / q16Scale := by
    apply Int.ediv_le_ediv hpos hmul
  rw [ofRawInt_toInt_eq_clamp, ofRawInt_toInt_eq_clamp]
  exact q16Clamp_monotone _ _ hdiv

/-- Addition is monotone in the left argument:
    if a ≤ b then a+c ≤ b+c. -/
theorem add_le_add (a b c : Q16_16) (h : a.toInt ≤ b.toInt) :
    (add a c).toInt ≤ (add b c).toInt := by
  unfold add
  simp [ofRawInt_toInt_eq_clamp]
  apply q16Clamp_monotone
  omega

/-- Absolute value of any Q16_16 value is non-negative. -/
theorem abs_nonneg (a : Q16_16) : (abs a).toInt ≥ 0 := by
  unfold abs neg
  have h := a.property.1
  split_ifs with hlt
  · rw [ofRawInt_toInt_eq_clamp]
    have h_nonneg : 0 ≤ -a.toInt := by omega
    exact q16Clamp_nonneg_of_nonneg h_nonneg
  · omega

/-- When addition does not saturate (the raw sum is in [q16MinRaw, q16MaxRaw]),
    the saturating `add` coincides with raw addition. Follows immediately
    from `q16Clamp_id_of_inRange` after rewriting via `ofRawInt_toInt_eq_clamp`. -/
theorem add_toInt_of_no_sat (a b : Q16_16)
    (hlo : q16MinRaw ≤ a.toInt + b.toInt) (hhi : a.toInt + b.toInt ≤ q16MaxRaw) :
    (add a b).toInt = a.toInt + b.toInt := by
  unfold add
  rw [ofRawInt_toInt_eq_clamp]
  exact q16Clamp_id_of_inRange _ hlo hhi

/-- When subtraction does not saturate, the saturating `sub` coincides with
    raw subtraction. -/
theorem sub_toInt_of_no_sat (a b : Q16_16)
    (hlo : q16MinRaw ≤ a.toInt - b.toInt) (hhi : a.toInt - b.toInt ≤ q16MaxRaw) :
    (sub a b).toInt = a.toInt - b.toInt := by
  unfold sub
  rw [ofRawInt_toInt_eq_clamp]
  exact q16Clamp_id_of_inRange _ hlo hhi

/-- Q16.16 multiplication rounds down to the integer floor for non-negative
    operands. The saturated result is at most `a*b/q16Scale`.

    The non-negative hypothesis gives `a*b/q16Scale ≥ 0 ≥ q16MinRaw`, so
    `q16Clamp` does not clamp the floor from below. Combined with the
    trivial `q16Clamp x ≤ q16MaxRaw` upper saturation, we get
    `q16Clamp (a*b/q16Scale) ≤ a*b/q16Scale`.

    Mirrors the `unfold mul; rw [ofRawInt_toInt_eq_clamp]` pattern used in
    `mul_mono_left`. -/
theorem mul_floor_le (a b : Q16_16) (ha : 0 ≤ a.toInt) (hb : 0 ≤ b.toInt) :
    (mul a b).toInt ≤ (a.toInt * b.toInt) / q16Scale := by
  unfold mul
  rw [ofRawInt_toInt_eq_clamp]
  have h_nonneg : 0 ≤ a.toInt * b.toInt / q16Scale := by
    apply Int.ediv_nonneg
    · nlinarith
    · norm_num [q16Scale]
  have h_ge_lo : q16MinRaw ≤ a.toInt * b.toInt / q16Scale := by
    dsimp [q16MinRaw]
    omega
  unfold q16Clamp
  by_cases h_hi : a.toInt * b.toInt / q16Scale > q16MaxRaw
  · dsimp [q16MaxRaw] at *
    simp [h_hi]
    omega
  · by_cases h_lo : a.toInt * b.toInt / q16Scale < q16MinRaw
    · dsimp [q16MinRaw] at *
      omega
    · simp [h_hi, h_lo]

/-- For non-negative operands, when the integer floor of the product fits
    within the Q16.16 range, multiplication is at least the floor.

    The hypothesis `a*b/q16Scale ≤ q16MaxRaw` rules out upper saturation,
    so `q16Clamp` does not clamp the floor from above. Together with
    `a*b/q16Scale ≥ 0 ≥ q16MinRaw` (from non-negativity) we get
    `q16Clamp (a*b/q16Scale) ≥ a*b/q16Scale`. -/
theorem mul_floor_ge (a b : Q16_16) (ha : 0 ≤ a.toInt) (hb : 0 ≤ b.toInt)
    (hhi : a.toInt * b.toInt / q16Scale ≤ q16MaxRaw) :
    (mul a b).toInt ≥ (a.toInt * b.toInt) / q16Scale := by
  unfold mul
  rw [ofRawInt_toInt_eq_clamp]
  have h_nonneg : 0 ≤ a.toInt * b.toInt / q16Scale := by
    apply Int.ediv_nonneg
    · nlinarith
    · norm_num [q16Scale]
  unfold q16Clamp
  by_cases h_hi : a.toInt * b.toInt / q16Scale > q16MaxRaw
  · dsimp [q16MaxRaw] at *
    omega
  · by_cases h_lo : a.toInt * b.toInt / q16Scale < q16MinRaw
    · dsimp [q16MinRaw] at *
      simp [h_hi, h_lo]
      omega
    · simp [h_hi, h_lo]

/-- Combined error bound for Q16.16 multiplication on non-negative operands.
    When the integer floor of the product fits in [q16MinRaw, q16MaxRaw],
    multiplication is exact: `(mul a b).toInt = a*b/q16Scale`.

    Combines `mul_floor_le` (upper bound) and `mul_floor_ge` (lower bound)
    via `Int.le_antisymm`. The two-sided bound closes the saturation
    envelope: no clamp from above (hhi) and no clamp from below
    (automatic from non-negativity). -/
theorem mul_floor_error (a b : Q16_16) (ha : 0 ≤ a.toInt) (hb : 0 ≤ b.toInt)
    (hhi : a.toInt * b.toInt / q16Scale ≤ q16MaxRaw) :
    (mul a b).toInt = (a.toInt * b.toInt) / q16Scale := by
  apply Int.le_antisymm
  · exact mul_floor_le a b ha hb
  · exact mul_floor_ge a b ha hb hhi

-- REMOVED: Theorems abs_mul_le and abs_triangle are FALSE for Q16_16 with floor division
-- abs_mul_le: counterexample a=3, b=-1 gives LHS=1, RHS=0
-- abs_triangle: counterexample a=3, b=-3 gives LHS=1, RHS=0
-- See TODO comments in original file for restructure guidance

/-- `ofNat` is monotone: a ≤ b → ofNat a ≤ ofNat b. -/
theorem ofNat_le (a b : Nat) (h : a ≤ b) : ofNat a ≤ ofNat b := by
  have h' : (a : Int) * q16Scale ≤ (b : Int) * q16Scale := by
    have h_nonneg : 0 ≤ (q16Scale : Int) := by norm_num [q16Scale]
    have h_int : (a : Int) ≤ (b : Int) := by exact_mod_cast h
    exact mul_le_mul_of_nonneg_right h_int h_nonneg
  exact ofRawInt_monotone _ _ h'

/-- `ofNat` returns non-negative values. -/
theorem ofNat_nonneg (n : Nat) : 0 ≤ (ofNat n).toInt := by
  unfold ofNat
  apply ofRawInt_toInt_nonneg
  have hpos : 0 ≤ (n : Int) * q16Scale :=
    mul_nonneg (Nat.cast_nonneg _) (by norm_num [q16Scale])
  exact hpos

/-- Addition is monotone in both arguments. -/
theorem add_le_add' (a b c d : Q16_16) (hac : a ≤ c) (hbd : b ≤ d) : a + b ≤ c + d := by
  have h_sum : a.toInt + b.toInt ≤ c.toInt + d.toInt := Int.add_le_add hac hbd
  have h_clamp : q16Clamp (a.toInt + b.toInt) ≤ q16Clamp (c.toInt + d.toInt) :=
    q16Clamp_monotone _ _ h_sum
  have h_add_toInt (x y : Q16_16) : (x + y).toInt = q16Clamp (x.toInt + y.toInt) := by
    rw [show x + y = add x y by rfl]
    rw [add]
    apply ofRawInt_toInt_eq_clamp
  calc
    (a + b).toInt = q16Clamp (a.toInt + b.toInt) := h_add_toInt _ _
    _ ≤ q16Clamp (c.toInt + d.toInt) := h_clamp
    _ = (c + d).toInt := (h_add_toInt _ _).symm

end Q16_16

-- ═══════════════════════════════════════════════════════════════════════════
-- Q0.64 signed normalized fraction
-- ═══════════════════════════════════════════════════════════════════════════

def q0_64MinRaw : Int := -9223372036854775808
def q0_64MaxRaw : Int := 9223372036854775807
def q0_64ScaleNat : Nat := 9223372036854775808

def q0_64ScaleFloat : Float := 9223372036854775808.0

/--
Q0.64 pure fraction representation.
The canonical proof model stores the signed raw integer in the Int64 range.
-/
abbrev Q0_64 := { x : Int // q0_64MinRaw ≤ x ∧ x ≤ q0_64MaxRaw }

instance : Repr Q0_64 where
  reprPrec q _ := repr q.val

instance : BEq Q0_64 where
  beq a b := a.val == b.val

instance : Inhabited Q0_64 where
  default := ⟨0, by constructor <;> norm_num [q0_64MinRaw, q0_64MaxRaw]⟩

instance : ToJson Q0_64 where
  toJson q := Json.mkObj [("val", toJson q.val)]

namespace Q0_64

@[ext]
theorem ext {a b : Q0_64} (h : a.val = b.val) : a = b := Subtype.ext h

@[inline]
def toInt (q : Q0_64) : Int := q.val

@[inline]
def ofRawInt (raw : Int) : Q0_64 :=
  if hhi : raw > q0_64MaxRaw then
    ⟨q0_64MaxRaw, by constructor <;> norm_num [q0_64MinRaw, q0_64MaxRaw]⟩
  else if hlo : raw < q0_64MinRaw then
    ⟨q0_64MinRaw, by constructor <;> norm_num [q0_64MinRaw, q0_64MaxRaw]⟩
  else
    ⟨raw, by
      constructor
      · dsimp [q0_64MinRaw, q0_64MaxRaw] at *; omega
      · dsimp [q0_64MinRaw, q0_64MaxRaw] at *; omega⟩

instance : FromJson Q0_64 where
  fromJson? j := do
    let raw : Int ← fromJson? (← j.getObjVal? "val")
    pure (ofRawInt raw)

/-- Maximum positive value. -/
def one : Q0_64 := ofRawInt q0_64MaxRaw

def zero : Q0_64 := ofRawInt 0

def ofRatio (num : Nat) (den : Nat) : Q0_64 :=
  if den = 0 then zero
  else ofRawInt (Int.ofNat (num * q0_64ScaleNat / den))

def half : Q0_64 := ofRawInt 4611686018427387903

def neg (x : Q0_64) : Q0_64 := ofRawInt (-x.toInt)
def add (a b : Q0_64) : Q0_64 := ofRawInt (a.toInt + b.toInt)
def sub (a b : Q0_64) : Q0_64 := ofRawInt (a.toInt - b.toInt)
def mul (a b : Q0_64) : Q0_64 :=
  ofRawInt ((a.toInt * b.toInt) / Int.ofNat q0_64ScaleNat)
def div (a b : Q0_64) : Q0_64 :=
  if b.toInt = 0 then one
  else ofRawInt ((a.toInt * Int.ofNat q0_64ScaleNat) / b.toInt)
def abs (x : Q0_64) : Q0_64 := if x.toInt < 0 then neg x else x

def ofFloat (f : Float) : Q0_64 :=
  if f.isNaN || f ≥ 1.0 then one
  else if f ≤ -1.0 then ofRawInt q0_64MinRaw
  else if f < 0.0 then
    ofRawInt (-(Int.ofNat ((-f * q0_64ScaleFloat).floor.toUInt64.toNat)))
  else
    ofRawInt (Int.ofNat ((f * q0_64ScaleFloat).floor.toUInt64.toNat))

def toFloat (q : Q0_64) : Float :=
  Float.ofInt q.toInt / q0_64ScaleFloat

instance : Add Q0_64 := ⟨add⟩
instance : Sub Q0_64 := ⟨sub⟩
instance : Mul Q0_64 := ⟨mul⟩
instance : Div Q0_64 := ⟨div⟩
instance : Neg Q0_64 := ⟨neg⟩

instance : LE Q0_64 where
  le a b := a.toInt ≤ b.toInt

instance : LT Q0_64 where
  lt a b := a.toInt < b.toInt

instance : DecidableRel (fun a b : Q0_64 => a ≤ b) :=
  fun a b => inferInstanceAs (Decidable (a.toInt ≤ b.toInt))

instance : DecidableRel (fun a b : Q0_64 => a < b) :=
  fun a b => inferInstanceAs (Decidable (a.toInt < b.toInt))

end Q0_64

-- ═══════════════════════════════════════════════════════════════════════════
-- ═══════════════════════════════════════════════════════════════════════════
-- Inverse Trig Witnesses
-- ═══════════════════════════════════════════════════════════════════════════

-- atan(0) = 0
#eval (Q16_16.atan Q16_16.zero).toInt  -- expect: 0

-- atan(1) = π/4 ≈ 0.7854 → raw ≈ 51471
#eval (Q16_16.atan Q16_16.one).toInt  -- expect: ~51471

-- atan(-1) = -π/4
#eval (Q16_16.atan (Q16_16.neg Q16_16.one)).toInt  -- expect: ~-51471

-- atan(large) ≈ π/2 (100.0 in Q16.16)
#eval (Q16_16.atan (Q16_16.ofRawInt 6553600)).toInt  -- expect: ~102943

-- asin(0) = 0
#eval (Q16_16.asin Q16_16.zero).toInt  -- expect: 0

-- asin(1) = π/2
#eval (Q16_16.asin Q16_16.one).toInt  -- expect: ~102943

-- asin(-1) = -π/2
#eval (Q16_16.asin (Q16_16.neg Q16_16.one)).toInt  -- expect: ~-102943

-- asin(0.5) ≈ 0.5236
#eval (Q16_16.asin (Q16_16.div Q16_16.one Q16_16.two)).toInt  -- expect: ~34306

-- acos(0) = π/2
#eval (Q16_16.acos Q16_16.zero).toInt  -- expect: ~102943

-- acos(1) = 0
#eval (Q16_16.acos Q16_16.one).toInt  -- expect: ~0

-- acos(-1) = π
#eval (Q16_16.acos (Q16_16.neg Q16_16.one)).toInt  -- expect: ~205887

-- atan2(1, 0) = π/2
#eval (Q16_16.atan2 Q16_16.one Q16_16.zero).toInt  -- expect: ~102943

-- atan2(0, 1) = 0
#eval (Q16_16.atan2 Q16_16.zero Q16_16.one).toInt  -- expect: 0

-- atan2(1, 1) = π/4
#eval (Q16_16.atan2 Q16_16.one Q16_16.one).toInt  -- expect: ~51471

-- atan2(-1, -1) = -3π/4
#eval (Q16_16.atan2 (Q16_16.neg Q16_16.one) (Q16_16.neg Q16_16.one)).toInt  -- expect: ~-154415

-- Pandigital π Approximation
-- ═══════════════════════════════════════════════════════════════════════════

namespace PandigitalPi

/-- High term: 3.8415926. -/
def highTerm : Q16_16 := Q16_16.ofRawInt 251819

/-- Low term: 0.7. -/
def lowTerm : Q16_16 := Q16_16.ofRawInt 45875

def piPandigital : Q16_16 := highTerm - lowTerm

def piDirect : Q16_16 := Q16_16.ofRawInt 205944

theorem piPandigitalCorrect : (piPandigital.toInt - piDirect.toInt).natAbs ≤ 1 := by
  native_decide

def spaceAnalysis : String :=
  "Pandigital pi: 6 bytes packed vs 4 bytes direct Q16.16 (trade-off for mathematical elegance)"

#eval piPandigital.toFloat
#eval piDirect.toFloat
#eval (piPandigital.toInt - piDirect.toInt).natAbs

end PandigitalPi

end SilverSight.FixedPoint

namespace SilverSight
  export FixedPoint (Q0_16 Q16_16 Q0_64)
  namespace Q16_16
    export FixedPoint.Q16_16
      (zero one negOne epsilon two infinity maxVal minVal ofNat satFromNat ofRatio toInt
       ofRawInt ofBits toBits ofFloat toFloat scale ofInt add sub mul div abs neg sqrt ln log2 pow sin log
       expNeg sat01 max min le ge gt lt recip ofRaw clip isNeg zero_mul mul_zero one_mul
       mul_one zero_add add_zero sub_self zero_toInt one_toInt epsilon_toInt
       epsilon_toInt_pos toInt_eq_zero_iff epsilon_add_pos zero_div mul_self_nonneg
       mul_toInt_nonneg ofRaw_toInt_nonneg mk_lt_half_nonneg add_one_omega_ge_one
       toInt_nonneg_le_maxVal add_pos_of_pos)
  end Q16_16
  namespace Q0_16
    export FixedPoint.Q0_16 (zero one half neg add sub mul div abs lt le gt ge toFloat ofFloat log2 min)
  end Q0_16
  namespace Q0_64
    export FixedPoint.Q0_64 (one zero ofRatio half neg add sub mul div abs toInt ofFloat toFloat)
  end Q0_64
  namespace PandigitalPi
    export FixedPoint.PandigitalPi (highTerm lowTerm piPandigital piDirect piPandigitalCorrect spaceAnalysis)
  end PandigitalPi
end SilverSight