/-
Copyright (c) 2026 Research Stack Contributors. All rights reserved.
Released under Apache 2.0 license.
-/
import SilverSight.FixedPoint
import Mathlib.Data.Finset.Basic
import Mathlib.NumberTheory.Divisors

/-! # Polynomial Factor Identity — Short-Sleeve Detection for RRC

This module ports Semantics.RRC.PolyFactorIdentity from the Research Stack
into SilverSight. It applies integer-to-polynomial decomposition as a
classification feature for the Rainbow Raccoon Compiler identity step.

The E8Sidon integration from the original file is replaced here by local,
minimal definitions of σ₃, σ₇, and the convolution LHS. The heavy E8/Sidon
theory is intentionally not pulled in; this file stays focused on the
poly-signature gate.

## Core Insight

When mathematical objects are accessed via zerocopy (mmap, shared memory,
framebuffer DMA), their raw limb structure is already exposed. Checking for
structured zero-blocks ("short-sleeve" patterns) is essentially free at that
boundary. When sparsity is detected, it flags the object for deeper
polynomial-based structural analysis.

## Integration

The `polySignature` function produces a `PolySignature` that RRC can use
alongside existing features (shape, alignment, receipt density) to classify
mathematical objects by their algebraic decomposability.

- High sparsity + low degree → likely factorizable (structured, "simple")
- Low sparsity + high degree → likely irreducible (complex, dense)
- High GCD of coefficients → common factor extractable (reducible)
-/

namespace SilverSight.RRC.PolyFactorIdentity

open SilverSight.FixedPoint
open SilverSight.FixedPoint.Q16_16
open Finset

-- ═══════════════════════════════════════════════════════════════════════════════
-- §0. Minimal E8 divisor-sum stubs (replaces the heavy Semantics.E8Sidon import)
-- ═══════════════════════════════════════════════════════════════════════════════

/-- Sum of divisors raised to the k-th power. -/
def sigmaK (k n : Nat) : Nat :=
  (Nat.divisors n).sum (fun d => d ^ k)

def sigma3 (n : Nat) : Nat := sigmaK 3 n
def sigma7 (n : Nat) : Nat := sigmaK 7 n

def convolutionLHS (n : Nat) : Nat :=
  (Finset.range (n - 1)).sum (fun j => sigma3 (j + 1) * sigma3 (n - j - 1))

-- ═══════════════════════════════════════════════════════════════════════════════
-- §1. Limb Decomposition (the zerocopy view)
-- ═══════════════════════════════════════════════════════════════════════════════

/-- Decompose a natural number into base-B limbs (least-significant first).
    This is the polynomial coefficient array: n = Σᵢ limbs[i] · Bⁱ. -/
def limbDecompose (n : Nat) (base : Nat) : List Nat :=
  if _h1 : base ≤ 1 then [n]
  else if _h2 : n = 0 then [0]
  else n % base :: limbDecompose (n / base) base
termination_by n
decreasing_by exact Nat.div_lt_self (Nat.pos_of_ne_zero _h2) (Nat.lt_of_not_le _h1)

-- Witnesses: limbs are LSB-first
#eval limbDecompose 2044 256    -- expect: [252, 7, 0] (7*256 + 252 = 2044)
#eval limbDecompose 65536 256   -- expect: [0, 0, 1, 0]
#eval limbDecompose 9 4         -- expect: [1, 2, 0] (2*4 + 1 = 9)

-- ═══════════════════════════════════════════════════════════════════════════════
-- §2. Coefficient Sparsity (the "short-sleeve" detector)
-- ═══════════════════════════════════════════════════════════════════════════════

/-- Count zero coefficients in a limb decomposition. -/
def zeroLimbCount (limbs : List Nat) : Nat :=
  limbs.filter (· == 0) |>.length

/-- Coefficient sparsity as Q16_16: (zeroCount / totalLimbs).
    Returns 0 if limbs is empty. -/
def coeffSparsity (limbs : List Nat) : Q16_16 :=
  if limbs.length == 0 then Q16_16.zero
  else Q16_16.ofRatio (zeroLimbCount limbs) limbs.length

/-- Short-sleeve detection threshold: 30% zero limbs triggers the flag.
    In Q16_16: 0.30 * 65536 ≈ 19661. -/
def shortSleeveThreshold : Q16_16 := Q16_16.ofRawInt 19661

/-- Does this integer exhibit short-sleeve structure when viewed in base-B? -/
def shortSleeveDetected (n : Nat) (base : Nat) : Bool :=
  let limbs := limbDecompose n base
  if limbs.length < 3 then false
  else (coeffSparsity limbs).toInt ≥ shortSleeveThreshold.toInt

-- Witnesses
#eval shortSleeveDetected 65536 256     -- expect: true (limbs [0,0,1] → 2/3 sparsity)
#eval shortSleeveDetected 2044 256      -- expect: false (limbs [252,7] → only 2 limbs)
#eval shortSleeveDetected 16777216 256  -- expect: true (limbs [0,0,0,1] → 3/4 sparsity)

-- ═══════════════════════════════════════════════════════════════════════════════
-- §3. Polynomial Signature (full structural features when flagged)
-- ═══════════════════════════════════════════════════════════════════════════════

/-- GCD of a list of natural numbers. -/
def listGcd : List Nat → Nat
  | [] => 0
  | [x] => x
  | x :: xs => Nat.gcd x (listGcd xs)

/-- Maximum coefficient in a limb decomposition. -/
def maxCoeff (limbs : List Nat) : Nat :=
  limbs.foldl Nat.max 0

/-- Sum of squares of coefficients — a proxy for "energy" of the polynomial. -/
def coeffEnergy (limbs : List Nat) : Nat :=
  limbs.foldl (fun acc c => acc + c * c) 0

/-- The polynomial degree: index of highest non-zero coefficient. -/
def polyDegree (limbs : List Nat) : Nat :=
  match limbs.reverse.dropWhile (· == 0) with
  | [] => 0
  | xs => xs.length - 1

/-- Full polynomial signature — the structural fingerprint. -/
structure PolySignature where
  base          : Nat
  degree        : Nat
  numTerms      : Nat
  sparsity      : Q16_16
  maxCoeffVal   : Nat
  gcdCoeffs     : Nat
  energy        : Nat
  isShortSleeve : Bool
  deriving Repr

/-- Compute the full polynomial signature for a natural number in base B. -/
def polySignature (n : Nat) (base : Nat) : PolySignature :=
  let limbs := limbDecompose n base
  let nonZero := limbs.filter (· != 0)
  { base          := base
    degree        := polyDegree limbs
    numTerms      := nonZero.length
    sparsity      := coeffSparsity limbs
    maxCoeffVal   := maxCoeff limbs
    gcdCoeffs     := listGcd nonZero
    energy        := coeffEnergy limbs
    isShortSleeve := shortSleeveDetected n base }

-- Witnesses
#eval polySignature 65536 256
#eval polySignature 16777216 256

-- ═══════════════════════════════════════════════════════════════════════════════
-- §4. E8Sidon Integration — Divisor Sum Polynomial Signatures
-- ═══════════════════════════════════════════════════════════════════════════════

def sigma3PolySig (n : Nat) (base : Nat) : PolySignature :=
  polySignature (sigma3 n) base

def sigma7PolySig (n : Nat) (base : Nat) : PolySignature :=
  polySignature (sigma7 n) base

def convPolySig (n : Nat) (base : Nat) : PolySignature :=
  polySignature (convolutionLHS n) base

#eval sigma7 4
#eval sigma7PolySig 4 256
#eval sigma7 6
#eval sigma7PolySig 6 256
#eval convolutionLHS 6
#eval convPolySig 6 256

-- ═══════════════════════════════════════════════════════════════════════════════
-- §5. RRC Feature Integration
-- ═══════════════════════════════════════════════════════════════════════════════

/-- Classification feature: polynomial decomposability score in [0, 1]. -/
def polyDecomposabilityScore (sig : PolySignature) : Q16_16 :=
  let sparsityComponent := (sig.sparsity.toInt * 26214) / q16Scale  -- * 0.4
  let gcdComponent : Int := if sig.gcdCoeffs > 1 then 19661 else 0  -- 0.3
  let termEfficiency : Int :=
    if sig.degree == 0 then 0
    else
      let ratio := (sig.numTerms * q16Scale) / (sig.degree + 1)
      let complement := q16Scale - ratio
      (complement * 19661) / q16Scale
  Q16_16.ofRawInt (sparsityComponent + gcdComponent + termEfficiency)

-- Witnesses
#eval polyDecomposabilityScore (polySignature 65536 256)
#eval polyDecomposabilityScore (polySignature 2044 256)
#eval polyDecomposabilityScore (polySignature 16777216 256)

-- ═══════════════════════════════════════════════════════════════════════════════
-- §6. Batch Scanner (zerocopy boundary integration)
-- ═══════════════════════════════════════════════════════════════════════════════

structure ZeroCopyScanResult where
  totalScanned    : Nat
  flaggedCount    : Nat
  flaggedIndices  : List Nat
  signatures      : List PolySignature
  deriving Repr

def zeroCopyScan (values : List Nat) (base : Nat) : ZeroCopyScanResult :=
  let indexed := values.zipIdx
  let flagged := indexed.filter (fun (v, _) => shortSleeveDetected v base)
  { totalScanned   := values.length
    flaggedCount   := flagged.length
    flaggedIndices := flagged.map Prod.snd
    signatures     := flagged.map (fun (v, _) => polySignature v base) }

#eval zeroCopyScan [2044, 65536, 255, 16777216, 42, 256] 256

-- ═══════════════════════════════════════════════════════════════════════════════
-- §7. Theorems — Polynomial Evaluation Correctness
-- ═══════════════════════════════════════════════════════════════════════════════

def polyEval : List Nat → Nat → Nat
  | [],         _    => 0
  | c :: rest,  base => c + base * polyEval rest base

/-- limbDecompose followed by polyEval recovers the original value. -/
theorem limbDecompose_polyEval_roundtrip (n : Nat) (base : Nat) (hb : base ≥ 2) :
    polyEval (limbDecompose n base) base = n := by
  induction n using Nat.strong_induction_on
  next n ih =>
  by_cases hn0 : n = 0
  · subst hn0
    unfold limbDecompose
    rw [dif_neg (show ¬(base ≤ 1) from by omega), dif_pos rfl]
    rfl
  · have h_div_lt : n / base < n := Nat.div_lt_self (by omega) (by omega)
    have hstep : limbDecompose n base = n % base :: limbDecompose (n / base) base := by
      conv_lhs => unfold limbDecompose; rw [dif_neg (show ¬(base ≤ 1) from by omega), dif_neg hn0]
    have hcons : polyEval (n % base :: limbDecompose (n / base) base) base =
                 n % base + base * polyEval (limbDecompose (n / base) base) base := rfl
    rw [hstep, hcons, ih (n / base) h_div_lt]
    exact Nat.mod_add_div n base

/-- Zero limbs bound the number of non-zero terms. -/
theorem zeroLimbs_bound_terms (limbs : List Nat) :
    (limbs.filter (· != 0)).length + zeroLimbCount limbs = limbs.length := by
  unfold zeroLimbCount
  induction limbs with
  | nil => simp
  | cons h t ih =>
    by_cases hh : h = 0
    · subst hh; simp; omega
    · simp [hh]; omega

/-- Cross-multiplication inequality for integer division. -/
lemma div_mul_div_le (a b c d : ℕ) (hb : b > 0) (hd : d > 0) (h : a * d ≥ c * b) : a / b ≥ c / d := by
  set k := c / d with hk
  have hc : k * d ≤ c := Nat.div_mul_le_self c d
  have ha_mul' : (k * d) * b ≤ a * d := by
    calc
      (k * d) * b ≤ c * b := Nat.mul_le_mul hc (le_refl _)
      _ ≤ a * d := h
  have ha : k * b ≤ a := by
    have htmp : (k * b) * d ≤ a * d := by
      calc
        (k * b) * d = (k * d) * b := by ring
        _ ≤ a * d := ha_mul'
    exact Nat.le_of_mul_le_mul_right htmp hd
  calc
    a / b ≥ (k * b) / b := Nat.div_le_div_right ha
    _ = k := by simp [hb.ne']
    _ = c / d := rfl

/-- Sparsity increases when a zero limb is prepended. -/
lemma sparsity_mono (z t : ℕ) (hz : z ≤ t) (h : z * 65536 / t ≥ 19661) : (z + 1) * 65536 / (t + 1) ≥ 19661 := by
  have ht0 : t > 0 := by
    by_contra! ht0
    have hz0 : t = 0 := by omega
    have : z * 65536 / t = 0 := by simp [hz0]
    rw [this] at h; omega
  have h_cross : (z + 1) * t ≥ z * (t + 1) := by
    calc
      (z + 1) * t = z * t + t := by ring
      _ ≥ z * t + z := Nat.add_le_add_left hz (z * t)
      _ = z * (t + 1) := by ring
  have h_ineq : (z + 1) * 65536 * t ≥ z * 65536 * (t + 1) := by
    calc
      (z + 1) * 65536 * t = 65536 * ((z + 1) * t) := by ring
      _ ≥ 65536 * (z * (t + 1)) := Nat.mul_le_mul_left _ h_cross
      _ = z * 65536 * (t + 1) := by ring
  have h_div : (z + 1) * 65536 / (t + 1) ≥ z * 65536 / t :=
    div_mul_div_le ((z + 1) * 65536) (t + 1) (z * 65536) t (by omega) ht0 h_ineq
  calc
    (z + 1) * 65536 / (t + 1) ≥ z * 65536 / t := h_div
    _ ≥ 19661 := h

lemma limbDecompose_mul_base (n : ℕ) (base : ℕ) (hb : base ≥ 2) (hn0 : n ≠ 0) :
    limbDecompose (n * base) base = 0 :: limbDecompose n base := by
  have hbase_gt1 : 1 < base := by omega
  have hbase0 : base > 0 := by omega
  have h_mul_nz : n * base ≠ 0 := mul_ne_zero hn0 (by omega)
  have h_mod : (n * base) % base = 0 := by simp
  have h_div : (n * base) / base = n := by
    simpa using Nat.mul_div_left n hbase0
  calc
    limbDecompose (n * base) base
        = (n * base) % base :: limbDecompose ((n * base) / base) base := by
          rw [limbDecompose]
          simp [hbase_gt1, h_mul_nz]
    _ = (0 :: limbDecompose n base) := by rw [h_mod, h_div]

lemma zeroLimbCount_le_length (limbs : List Nat) : zeroLimbCount limbs ≤ limbs.length := by
  unfold zeroLimbCount
  exact List.length_filter_le (· == 0) limbs

lemma coeffSparsity_val (limbs : List Nat) (hlen : limbs.length > 0) :
    (coeffSparsity limbs).toInt = (zeroLimbCount limbs * 65536 / limbs.length : ℤ) := by
  unfold coeffSparsity
  have hnonzero : limbs.length ≠ 0 := by omega
  simp [hnonzero, Q16_16.ofRatio, ofRawInt_toInt_eq_clamp]
  have hz : zeroLimbCount limbs ≤ limbs.length := zeroLimbCount_le_length _
  have hmul : zeroLimbCount limbs * 65536 ≤ limbs.length * 65536 := Nat.mul_le_mul_right _ hz
  have hdiv : zeroLimbCount limbs * 65536 / limbs.length ≤ limbs.length * 65536 / limbs.length :=
    Nat.div_le_div_right hmul
  have hcancel : limbs.length * 65536 / limbs.length = 65536 := by
    simpa [Nat.mul_comm] using Nat.mul_div_cancel 65536 hlen
  rw [hcancel] at hdiv
  apply q16Clamp_id_of_inRange
  · have h_nonneg_z : (0 : ℤ) ≤ (zeroLimbCount limbs : ℤ) := Nat.cast_nonneg _
    have h_nonneg_65536 : (0 : ℤ) ≤ (65536 : ℤ) := by norm_num
    have h_nonneg_prod : (0 : ℤ) ≤ (zeroLimbCount limbs : ℤ) * (65536 : ℤ) :=
      mul_nonneg h_nonneg_z h_nonneg_65536
    have h_nonneg_div : (0 : ℤ) ≤ (zeroLimbCount limbs : ℤ) * (65536 : ℤ) / (limbs.length : ℤ) :=
      Int.ediv_nonneg h_nonneg_prod (Nat.cast_nonneg _)
    unfold q16MinRaw
    calc
      (-2147483648 : ℤ) ≤ (0 : ℤ) := by norm_num
      _ ≤ (zeroLimbCount limbs : ℤ) * (65536 : ℤ) / (limbs.length : ℤ) := h_nonneg_div
  · have hdiv_int : (zeroLimbCount limbs * 65536 / limbs.length : ℤ) ≤ (65536 : ℤ) := by
      exact_mod_cast hdiv
    have h_65536_max : (65536 : ℤ) ≤ q16MaxRaw := by unfold q16MaxRaw; norm_num
    exact le_trans hdiv_int h_65536_max

/-- shortSleeveDetected is monotone under multiplication by base. -/
theorem shortSleeve_mono_zero_prepend (n : Nat) (base : Nat) (hb : base ≥ 2)
    (h : shortSleeveDetected n base = true) :
    shortSleeveDetected (n * base) base = true := by
  simp only [shortSleeveDetected] at h
  split_ifs at h with hlt
  simp only [decide_eq_true_iff, ge_iff_le] at h
  push Not at hlt
  have hn0 : n ≠ 0 := by
    intro heq; subst heq
    have heq0 : limbDecompose 0 base = [0] := by
      conv_lhs => unfold limbDecompose; rw [dif_neg (by omega : ¬(base ≤ 1)), dif_pos rfl]
    simp [heq0] at hlt
  have hmul := limbDecompose_mul_base n base hb hn0
  set limbs := limbDecompose n base with hlimbs
  set z := zeroLimbCount limbs with hz_def
  set t := limbs.length with ht_def
  have hlen_pos : t > 0 := by omega
  have hcs  := coeffSparsity_val limbs hlen_pos
  have hth  : shortSleeveThreshold.toInt = 19661 := by rfl
  have hsp_nat : z * 65536 / t ≥ 19661 := by
    have hℤ : (19661 : ℤ) ≤ (z : ℤ) * 65536 / t := calc
      (19661 : ℤ) = shortSleeveThreshold.toInt := hth.symm
      _ ≤ (coeffSparsity limbs).toInt          := h
      _ = (z : ℤ) * 65536 / t                 := hcs
    exact_mod_cast hℤ
  have hz_le  := zeroLimbCount_le_length limbs
  have hmono  := sparsity_mono z t hz_le hsp_nat
  have hzero' : zeroLimbCount (0 :: limbs) = z + 1 := by simp [zeroLimbCount, hz_def]
  have hlen'  : (0 :: limbs).length = t + 1 := by simp [ht_def]
  have hcs_new := coeffSparsity_val (0 :: limbs) (by simp)
  rw [hzero', hlen'] at hcs_new
  simp only [shortSleeveDetected, hmul]
  split_ifs with hlt'
  · simp only [hlen'] at hlt'; omega
  · simp only [decide_eq_true_iff, ge_iff_le]
    rw [hcs_new, hth]
    exact_mod_cast hmono

end SilverSight.RRC.PolyFactorIdentity
