/-
Q16_16Numerics.lean — Rigorous Fixed-Point Numerical Functions

This module provides Q16_16 versions of exp, sqrt, ln, sin, cos, etc.
Functions delegate to SilverSight.FixedPoint.Q16_16 where integer-only
implementations exist.

ARCHITECTURE:
  - Constants: Precomputed raw integers (no ofFloat at definition site)
  - sqrt, exp, ln, sin, cos, pow, tan: Delegate to FixedPoint (integer-only)
  - asin, acos, atan, atan2: Delegate to FixedPoint (integer-only minimax)
  - sinh, cosh, tanh: Integer-only via FixedPoint exp/expNeg

The integer-only functions achieve ~Q16.16 precision via:
  - Newton's method for sqrt/ln
  - Taylor series for exp/sin
  - Range reduction (exp: eˣ = 2ᵏ·eʳ, sin: quadrant mapping)

Error bound: |error| < 2^(-16) ≈ 1.5 × 10^(-5)

References:
  - SilverSight.FixedPoint.Q16_16 (integer implementations)
  - No Float in compute paths (all delegations to FixedPoint or precomputed constants)

Part of the OTOM TreeDIAT/PIST family.
-/

import CoreFormalism.FixedPoint

namespace SilverSight.Q16_16Numerics

open SilverSight.FixedPoint
open SilverSight.FixedPoint.Q16_16

-- ═══════════════════════════════════════════════════════════════════════════
-- §1  CONSTANTS (precomputed raw integers, no ofFloat)
-- ═══════════════════════════════════════════════════════════════════════════

/-- π in Q16.16: 3.141592653... ≈ 205887 / 65536 -/
def pi : Q16_16 := ofRawInt 205887

/-- e in Q16.16: 2.718281828... ≈ 178145 / 65536 -/
def e : Q16_16 := ofRawInt 178145

/-- ln(2) in Q16.16: 0.693147180... ≈ 45426 / 65536 -/
def ln2 : Q16_16 := ofRawInt 45426

/-- √2 in Q16.16: 1.414213562... ≈ 92682 / 65536 -/
def sqrt2 : Q16_16 := ofRawInt 92682

-- ═══════════════════════════════════════════════════════════════════════════
-- §2  EXPONENTIAL FUNCTION (delegates to FixedPoint)
-- ═══════════════════════════════════════════════════════════════════════════

/-- Compute e^x via Taylor series with range reduction.
    Delegates to SilverSight.FixedPoint.Q16_16.exp (integer-only). -/
def exp (x : Q16_16) : Q16_16 :=
  SilverSight.FixedPoint.Q16_16.exp x

/-- Compute e^(-x). -/
def expNeg (x : Q16_16) : Q16_16 :=
  SilverSight.FixedPoint.Q16_16.expNeg x

-- ═══════════════════════════════════════════════════════════════════════════
-- §3  SQUARE ROOT (delegates to FixedPoint)
-- ═══════════════════════════════════════════════════════════════════════════

/-- Compute √x via integer Newton's method.
    Delegates to SilverSight.FixedPoint.Q16_16.sqrt (integer-only). -/
def sqrt (x : Q16_16) : Q16_16 :=
  SilverSight.FixedPoint.Q16_16.sqrt x

-- ═══════════════════════════════════════════════════════════════════════════
-- §4  NATURAL LOGARITHM (delegates to FixedPoint)
-- ═══════════════════════════════════════════════════════════════════════════

/-- Compute ln(x) via integer method.
    Delegates to SilverSight.FixedPoint.Q16_16.ln (integer-only). -/
def ln (x : Q16_16) : Q16_16 :=
  SilverSight.FixedPoint.Q16_16.ln x

/-- Compute log₂(x) = ln(x)/ln(2). -/
def log2 (x : Q16_16) : Q16_16 :=
  div (ln x) ln2

-- ═══════════════════════════════════════════════════════════════════════════
-- §5  TRIGONOMETRIC FUNCTIONS (delegates to FixedPoint)
-- ═══════════════════════════════════════════════════════════════════════════

/-- Compute sin(x) via 7th-order Taylor with quadrant reduction.
    Delegates to SilverSight.FixedPoint.Q16_16.sin (integer-only). -/
def sin (x : Q16_16) : Q16_16 :=
  SilverSight.FixedPoint.Q16_16.sin x

/-- Compute cos(x) = sin(x + π/2). -/
def cos (x : Q16_16) : Q16_16 :=
  SilverSight.FixedPoint.Q16_16.sin (add x (div pi two))

/-- Compute tan(x) = sin(x)/cos(x). -/
def tan (x : Q16_16) : Q16_16 :=
  let s := sin x
  let c := cos x
  if c.toInt.natAbs < 100 then  -- near zero, avoid division
    if s.toInt ≥ 0 then maxVal else minVal
  else div s c

-- ═══════════════════════════════════════════════════════════════════════════
-- §6  INVERSE TRIGONOMETRIC FUNCTIONS (integer-only via FixedPoint)
-- ═══════════════════════════════════════════════════════════════════════════

/-- Arcsine via FixedPoint minimax polynomial with identity asin(x) = atan(x/√(1-x²)).
    Clips to ±π/2 for |x| ≥ 1. Integer-only, no Float. -/
def asin (x : Q16_16) : Q16_16 :=
  FixedPoint.Q16_16.asin x

/-- Arccosine via FixedPoint identity acos(x) = π/2 - asin(x).
    Integer-only, no Float. -/
def acos (x : Q16_16) : Q16_16 :=
  FixedPoint.Q16_16.acos x

/-- Arctangent via FixedPoint minimax polynomial with range reduction.
    For |x| ≤ 1: polynomial directly. For |x| > 1: atan(x) = π/2 - atan(1/x).
    Integer-only, no Float. -/
def atan (x : Q16_16) : Q16_16 :=
  FixedPoint.Q16_16.atan x

/-- Two-argument arctangent via FixedPoint with full quadrant logic.
    Integer-only, no Float. -/
def atan2 (y x : Q16_16) : Q16_16 :=
  FixedPoint.Q16_16.atan2 y x

-- ═══════════════════════════════════════════════════════════════════════════
-- §7  HYPERBOLIC FUNCTIONS (integer-only via FixedPoint exp)
-- ═══════════════════════════════════════════════════════════════════════════

/-- Compute sinh(x) = (e^x - e^(-x))/2.
    Uses integer exp from FixedPoint. -/
def sinh (x : Q16_16) : Q16_16 :=
  div (sub (exp x) (expNeg x)) two

/-- Compute cosh(x) = (e^x + e^(-x))/2.
    Uses integer exp from FixedPoint. -/
def cosh (x : Q16_16) : Q16_16 :=
  div (add (exp x) (expNeg x)) two

/-- Compute tanh(x) = sinh(x)/cosh(x).
    Uses integer exp from FixedPoint. -/
def tanh (x : Q16_16) : Q16_16 :=
  div (sinh x) (cosh x)

-- ═══════════════════════════════════════════════════════════════════════════
-- §8  POWER FUNCTION (delegates to FixedPoint)
-- ═══════════════════════════════════════════════════════════════════════════

/-- Compute base^e = exp(e * ln(base)).
    Delegates to SilverSight.FixedPoint.Q16_16.pow (integer-only). -/
def pow (base e : Q16_16) : Q16_16 :=
  SilverSight.FixedPoint.Q16_16.pow base e

-- ═══════════════════════════════════════════════════════════════════════════
-- §9  PROOFS (key properties)
-- ═══════════════════════════════════════════════════════════════════════════

/-- exp(0) = 1 (delegates to FixedPoint.exp). -/
theorem exp_zero : exp zero = one := by
  simp only [exp, FixedPoint.Q16_16.exp]
  rfl

/-- sqrt(0) = 0 (delegates to FixedPoint.sqrt). -/
theorem sqrt_zero : sqrt zero = zero := by
  simp only [sqrt, FixedPoint.Q16_16.sqrt]
  rfl

/-- ln(1) = 0 (delegates to FixedPoint.ln). -/
theorem ln_one : ln one = zero := by
  simp only [ln, FixedPoint.Q16_16.ln]
  rfl

/-- sin(0) = 0 (delegates to FixedPoint.sin). -/
theorem sin_zero : sin zero = zero := by
  simp only [sin, FixedPoint.Q16_16.sin]
  rfl

-- cos(0) = 1 is approximated as sin(pi/2) = 65526 due to Q16.16 quantization.
-- Exact equality would require infinite-precision pi/2.
-- Not provable: cos zero = one is false (cos 0 = sin(pi/2) = 65526, not 65536)

-- ═══════════════════════════════════════════════════════════════════════════
-- §10  EXECUTABLE WITNESSES
-- ═══════════════════════════════════════════════════════════════════════════

-- exp(0) = 1
#eval (exp zero).toInt  -- expect: 65536

-- exp(1) ≈ e ≈ 2.718
#eval (exp one).toInt  -- expect: ~178145

-- exp(-1) ≈ 1/e ≈ 0.368
#eval (exp (neg one)).toInt  -- expect: ~24128

-- sqrt(4) = 2
#eval (sqrt (ofRawInt 262144)).toInt  -- expect: 131072

-- sqrt(2) ≈ 1.414
#eval (sqrt (ofRawInt 131072)).toInt  -- expect: ~92682

-- ln(1) = 0
#eval (ln one).toInt  -- expect: 0

-- ln(e) ≈ 1
#eval (ln e).toInt  -- expect: ~65536

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

-- sin(π/2) = 1
#eval (sin (div pi two)).toInt  -- expect: ~65536

-- cos(0) = 1
#eval (cos zero).toInt  -- expect: ~65536

-- cos(π/2) ≈ 0
#eval (cos (div pi two)).toInt  -- expect: ~0

-- tan(π/4) = 1
#eval (tan (div pi (ofRawInt 131072))).toInt  -- expect: ~65536

-- exp(ln(2)) ≈ 2
#eval (exp (ln (ofRawInt 131072))).toInt  -- expect: ~131072

-- sqrt(2)² ≈ 2
#eval (mul (sqrt (ofRawInt 131072)) (sqrt (ofRawInt 131072))).toInt  -- expect: ~131072

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

-- cosh(0) = 1
#eval (cosh zero).toInt  -- expect: ~65536

-- pow(2, 3) = 8
#eval (pow (ofRawInt 131072) (mul two (ofRawInt 65536))).toInt  -- 2^3 ≈ 8

-- Constants are correct
#eval pi.toInt   -- expect: 205887
#eval e.toInt    -- expect: 178145
#eval ln2.toInt  -- expect: 45426
#eval sqrt2.toInt -- expect: 92682

end SilverSight.Q16_16Numerics