SilverSight/julia/CoreFormalism/Q16_16.jl

"""
    SilverSight.FixedPoint.Q16_16 — Julia Port

Mirrors `Core/SilverSight/FixedPoint.lean` (Lean 4).

Design rules:
  * Semantic value is a signed bounded `Int32` raw integer.
  * Saturation is performed by `ofRawInt`.
  * Float is external boundary only (`ofFloat`, `toFloat`).

All operations match Lean semantics: saturation, rounding mode,
and constant values are verified against the Lean `#eval` witnesses
in `Q16_16Numerics.lean`.
"""
module Q16_16

export Q16_16Value, Q16_SCALE,
       of_raw_int, to_int,
       zero, one, epsilon, two, half,
       q16_min_raw, q16_max_raw, q16_scale,
       add, sub, mul, div_op, neg, abs_,
       sqrt, ln_nat, exp, exp_neg, sin, pow,
       of_float, to_float,
       lt, le, gt, ge,
       q16_pct1, q16_pct70, q16_pct30, q16_150pct

const q16_min_raw = -2147483648
const q16_max_raw = 2147483647
const q16_scale   = 65536
const Q16_SCALE   = q16_scale

"""
    Q16_16Value

Signed Q16.16 fixed-point raw integer, saturated into `[q16_min_raw, q16_max_raw]`.
Concrete type: `Int32` (matches `ofBits`/`toBits` hardware representation).
"""
const Q16_16Value = Int32

"""
    clamp_raw(raw::Int) -> Int

Saturating clamp matching Lean `q16Clamp`. Returns `q16_min_raw` if `raw < q16_min_raw`,
`q16_max_raw` if `raw > q16_max_raw`, else `raw`.
"""
function clamp_raw(raw::Int)::Int
    if raw > q16_max_raw
        return q16_max_raw
    elseif raw < q16_min_raw
        return q16_min_raw
    else
        return raw
    end
end

"""
    of_raw_int(raw) -> Q16_16Value

Match Lean `Q16_16.ofRawInt`: saturating clamp into valid range.
"""
function of_raw_int(raw::Integer)::Q16_16Value
    Q16_16Value(clamp_raw(Int(raw)))
end

of_raw_int(raw::Int32)::Q16_16Value = raw

"""
    to_int(q) -> Int

Match Lean `Q16_16.toInt`.
"""
to_int(q::Q16_16Value)::Int = Int(q)

# ── Constants (match Lean `#eval` witnesses exactly) ────────────────────────

const zero    = Q16_16Value(0)
const one     = Q16_16Value(65536)
const neg_one = Q16_16Value(-65536)
const epsilon = Q16_16Value(1)         # Lean: epsilon = ⟨1, ...⟩
const two     = Q16_16Value(131072)    # 2 * q16_scale
const half    = Q16_16Value(32768)     # q16_scale / 2
const max_val = Q16_16Value(q16_max_raw)
const min_val = Q16_16Value(q16_min_raw)
const infinity = max_val

# NUVMAP constants (from projection_engine.py)
const q16_pct1   = Q16_16Value(655)    # 0.01 floor
const q16_pct70  = Q16_16Value(45875)  # 0.7 floor
const q16_pct30  = Q16_16Value(19661)  # 0.3 ceil (banker's)
const q16_150pct = Q16_16Value(98304)  # 1.5 exact

# Mathemtical constants (from Q16_16Numerics.lean)
const pi_val    = Q16_16Value(205887)  # π
const e_val     = Q16_16Value(178145)  # e
const ln2_val   = Q16_16Value(45426)   # ln(2)
const sqrt2_val = Q16_16Value(92682)   # √2

# ── Arithmetic (saturating, matching Lean) ──────────────────────────────────

"""
    add(a, b)

Saturating addition. Matches Lean `Q16_16.add`.
"""
add(a::Q16_16Value, b::Q16_16Value)::Q16_16Value =
    of_raw_int(to_int(a) + to_int(b))

"""
    sub(a, b)

Saturating subtraction. Matches Lean `Q16_16.sub`.
"""
sub(a::Q16_16Value, b::Q16_16Value)::Q16_16Value =
    of_raw_int(to_int(a) - to_int(b))

"""
    mul(a, b)

Saturating Q16.16 multiplication with banker's rounding: `round((a * b) / 65536)`.
Matches Python `q_mul` semantics. Uses `Int64` intermediate.
"""
function mul(a::Q16_16Value, b::Q16_16Value)::Q16_16Value
    prod = Int64(to_int(a)) * Int64(to_int(b))
    of_raw_int(Int(round(prod / q16_scale)))
end

"""
    div_op(a, b)

Saturating Q16.16 division with banker's rounding: `round((a * 65536) / b)`.
Matches Python `q_div` semantics. Returns `infinity` if `b == 0`.
Uses `Int64` intermediate.
"""
function div_op(a::Q16_16Value, b::Q16_16Value)::Q16_16Value
    if to_int(b) == 0
        return infinity
    end
    prod = Int64(to_int(a)) * Int64(q16_scale)
    of_raw_int(Int(round(prod / Int64(to_int(b)))))
end

"""
    neg(a)

Saturating negation. Matches Lean `Q16_16.neg`.
"""
neg(a::Q16_16Value)::Q16_16Value = of_raw_int(-to_int(a))

"""
    abs_(a)

Absolute value. Matches Lean `Q16_16.abs`.
"""
abs_(a::Q16_16Value)::Q16_16Value = to_int(a) < 0 ? neg(a) : a

# ── Comparisons ─────────────────────────────────────────────────────────────

lt(a::Q16_16Value, b::Q16_16Value)::Bool = to_int(a) < to_int(b)
le(a::Q16_16Value, b::Q16_16Value)::Bool = to_int(a) ≤ to_int(b)
gt(a::Q16_16Value, b::Q16_16Value)::Bool = to_int(a) > to_int(b)
ge(a::Q16_16Value, b::Q16_16Value)::Bool = to_int(a) ≥ to_int(b)

# ── Float boundary (external only) ──────────────────────────────────────────

"""
    of_float(f)

Float → Q16_16. Matches Lean `Q16_16.ofFloat`:
  - NaN or ≥ 32768.0 → infinity
  - ≤ -32768.0 → min_val
  - floor(f * 65536) with sign handling

External boundary only.
"""
function of_float(f::Float64)::Q16_16Value
    if isnan(f) || f ≥ 32768.0
        return infinity
    elseif f ≤ -32768.0
        return min_val
    else
        return of_raw_int(Int(floor(f * 65536.0)))
    end
end

of_float(f::Float32)::Q16_16Value = of_float(Float64(f))

"""
    to_float(q)

Q16_16 → Float64. Matches Lean `Q16_16.toFloat`.
"""
to_float(q::Q16_16Value)::Float64 = Float64(to_int(q)) / 65536.0

# ── Numeric functions ───────────────────────────────────────────────────────

"""
    int_sqrt(n)

Integer square root via Newton's method (floor). Matches Lean `intSqrt`.
Fuel-limited to 64 iterations.
"""
function int_sqrt(n::Int)::Int
    n ≤ 0 && return 0
    x = div(n, 2) + 1
    for _ in 1:64
        x_new = div(x + div(n, x), 2)
        x_new ≥ x && return x
        x = x_new
    end
    return x
end

"""
    sqrt(q)

Q16.16 square root via integer Newton. Matches Lean `Q16_16.sqrt`.
"""
sqrt(q::Q16_16Value)::Q16_16Value =
    to_int(q) ≤ 0 ? zero : of_raw_int(int_sqrt(to_int(q) * q16_scale))

"""
    ln_nat(q)

Natural logarithm approximation around 1.0. Matches Lean `Q16_16.ln`.
"""
function ln_nat(q::Q16_16Value)::Q16_16Value
    x = to_int(q)
    x ≤ 0 && return zero
    y = x - q16_scale
    y2 = div(y * y, q16_scale)
    y3 = div(y * y2, q16_scale)
    of_raw_int(y - div(y2, 2) + div(y3, 3))
end

"""
    exp_neg(x)

e^(-x) piecewise approximation. Matches Lean `Q16_16.expNeg`.
"""
function exp_neg(x::Q16_16Value)::Q16_16Value
    rx = to_int(x)
    if rx ≥ 0x00030000; return zero
    elseif rx ≥ 0x00020000; return of_raw_int(0x00004D29)
    elseif rx ≥ 0x00010000; return of_raw_int(0x0000C5C0)
    else return of_raw_int(0x0001C5C0)
    end
end

"""
    exp(x)

Q16.16 exponential via Taylor with range reduction eˣ = 2ᵏ·eʳ.
Matches Lean `Q16_16.exp`.
"""
function exp(x::Q16_16Value)::Q16_16Value
    rx = to_int(x)
    if rx ≤ -q16_scale; return zero
    elseif rx ≥ 4 * q16_scale; return max_val
    end

    ln2_raw = 45426
    k = div(rx, ln2_raw)
    r = rx - k * ln2_raw
    r2 = div(r * r, q16_scale)
    r3 = div(r * r2, q16_scale)
    r4 = div(r2 * r2, q16_scale)
    r5 = div(r2 * r3, q16_scale)
    r6 = div(r3 * r3, q16_scale)
    taylor = q16_scale + r + div(r2, 2) + div(r3, 6) + div(r4, 24) + div(r5, 120) + div(r6, 720)

    shifted = if k ≥ 15; q16_max_raw
    elseif k ≤ -15; 0
    elseif k ≥ 0; taylor << UInt(k)
    else; taylor >> UInt(-k)
    end

    of_raw_int(shifted)
end

"""
    sin(x)

Q16.16 sine via 7th-order Taylor with quadrant reduction.
Matches Lean `Q16_16.sin`.
"""
function sin(x::Q16_16Value)::Q16_16Value
    pi_raw = 205887
    two_pi_raw = 411774
    raw = to_int(x) % two_pi_raw
    raw = raw < 0 ? raw + two_pi_raw : raw
    half_pi = div(pi_raw, 2)

    q = if raw < half_pi; 0
    elseif raw < pi_raw; 1
    elseif raw < pi_raw + half_pi; 2
    else; 3
    end

    reduced = if q == 0; raw
    elseif q == 1; pi_raw - raw
    elseif q == 2; raw - pi_raw
    else; two_pi_raw - raw
    end

    t = reduced
    t2 = div(t * t, q16_scale)
    t3 = div(t * t2, q16_scale)
    t5 = div(t3 * t2, q16_scale)
    t7 = div(t5 * t2, q16_scale)
    sin_pos = t - div(t3, 6) + div(t5, 120) - div(t7, 5040)
    result = q < 2 ? sin_pos : -sin_pos
    of_raw_int(max(q16_min_raw, min(q16_max_raw, result)))
end

"""
    pow(base, e)

Q16.16 power: base^e = exp(e * ln(base)). Matches Lean `Q16_16.pow`.
"""
function pow(base::Q16_16Value, e::Q16_16Value)::Q16_16Value
    to_int(base) ≤ 0 && return to_int(e) == 0 ? one : zero
    to_int(e) == 0 && return one
    to_int(e) == q16_scale && return base
    return exp(mul(ln_nat(base), e))
end

# ── NUVMAP convenience ──────────────────────────────────────────────────────

"""
    q_mul(a, b) — Q16.16 multiplication with banker's rounding.
    q_div(a, b) — Q16.16 division with banker's rounding. Returns 0 if b == 0.
    q_add(a, b) — Saturating addition.
    q_sub(a, b) — Saturating subtraction.

These match the naming in `python/nuvmap/projection_engine.py` and
`tests/test_nuvmap_equivalence.py`.
"""
q_mul(a::Q16_16Value, b::Q16_16Value)::Q16_16Value = mul(a, b)
q_div(a::Q16_16Value, b::Q16_16Value)::Q16_16Value = div_op(a, b)
q_add(a::Q16_16Value, b::Q16_16Value)::Q16_16Value = add(a, b)
q_sub(a::Q16_16Value, b::Q16_16Value)::Q16_16Value = sub(a, b)

# ── Instance methods for convenience ──────────────────────────────────────

Base.:+(a::Q16_16Value) = a
Base.:-(a::Q16_16Value) = neg(a)
Base.:+(a::Q16_16Value, b::Q16_16Value) = add(a, b)
Base.:-(a::Q16_16Value, b::Q16_16Value) = sub(a, b)
Base.:*(a::Q16_16Value, b::Q16_16Value) = mul(a, b)
Base.:/(a::Q16_16Value, b::Q16_16Value) = div_op(a, b)
Base.:<(a::Q16_16Value, b::Q16_16Value) = lt(a, b)
Base.:<=(a::Q16_16Value, b::Q16_16Value) = le(a, b)
Base.:>(a::Q16_16Value, b::Q16_16Value) = gt(a, b)
Base.:>=(a::Q16_16Value, b::Q16_16Value) = ge(a, b)
Base.abs(a::Q16_16Value) = abs_(a)
Base.sqrt(a::Q16_16Value) = sqrt(a)
Base.exp(a::Q16_16Value) = exp(a)
Base.sin(a::Q16_16Value) = sin(a)
Base.:^(a::Q16_16Value, b::Q16_16Value) = pow(a, b)

end