Upload 6 files

q3_k.py

@@ -1,22 +1,41 @@
-"""Full HPC V2 Q3_K quantizer.
-
-Implements the HPC-Quantize V2 decision stack for llama.cpp/GGML ``Q3_K``
-(256 weights per 110-byte block = 3.4375 bpw):
-
-1. **Geometry-Preserving Analog + WLS Seeds**:
-   Finds an optimal signed scale from the raw weights and an iso-moment
-   companded analog for every 16-weight sub-block.
-2. **24-Candidate Scale Lattice & Extended Spectral Objective**:
-   Evaluates 24 scale multipliers (mapped to 6 symbolic quhit states) on
-   importance-weighted SSE + D6/Vesica fold pyramid + DC penalty, including
-   negative-scale sign-flip candidates.
-3. **6-State Sequential Sieve & 6-State Viterbi DP**:
-   Runs the 6-state sequential Sieve with parity back-action and a 6-state
-   Viterbi dynamic program across the sub-block chain.
-4. **Asymmetric Sign-Flip fp16-ULP Polish & Closed-Form Vesica/DC Shaping**:
-   Polishes each sub-block scale over ``+/-8`` fp16 ULPs, then jointly fits the
-   GGML super-scale ``d`` and signed 6-bit sub-block scales ``sc``.  The final
-   3-bit codes ``[-4, +3]`` are emitted in the exact ``block_q3_K`` layout.
+"""Full HPC V2 Q3_K quantizer (CUDA GPU + Vectorized CPU).
+
+Implements the complete HPC-Quantize V2 decision stack for llama.cpp/GGML
+``Q3_K`` (256 weights per 110-byte block = 3.4375 bpw).  The Q3_K hierarchy is
+the Q4_0 pipeline applied twice:
+
+* **Level A — 16-weight sub-blocks** (the per-sub-block scale ``s_j`` search):
+  runs the exact Q4_0 decision stack on every 16-weight sub-block with the
+  signed 3-bit grid ``[-4, +3]``:
+  1. Geometry-Preserving Analog + WLS seeds (``+d`` and asymmetric ``-d``).
+  2. 24-candidate scale lattice & extended spectral objective.
+  3. 6-state sequential Sieve & 6-state Viterbi DP across the sub-block chain.
+  4. Asymmetric sign-flip fp16-ULP polish.
+  5. Closed-form 1-D quadratic Vesica + DC code projection (``sse_budget``).
+
+* **Level B — 256-weight super-blocks** (the GGML ``(d, sc)`` fit, no
+  shortcuts relative to Q4_0):
+  1. **Geometry-Preserving Analog + WLS Seeds** for the super-scale ``d``:
+     four seeds (``+max``/``-max`` grid midpoint ``/31.5`` and the iso-moment
+     companded analog ``+/-``), each refined by 4 rounds of weighted least
+     squares over the 6-bit sub-block scales with per-sub-block importance,
+     fp16-rounded and scored by the true 3-bit encode SSE (``+d`` and ``-d``).
+  2. **24-Candidate Super-Scale Lattice & Extended Spectral Objective**:
+     ``cand_d = fp16(wls_d * MULTS_D)``, each candidate jointlyLS derives the
+     signed 6-bit scales ``sc ∈ [-32, 31]`` and 3-bit codes ``[-4, +3]`` and is
+     scored on importance-weighted SSE + D6/Vesica fold pyramid + DC penalty
+     over the full 256-weight block.
+  3. **6-State Sequential Sieve & 6-State Viterbi DP** over the 24-candidate
+     lattice with free-energy bin reduction and greedy override.
+  4. **Asymmetric Sign-Flip fp16-ULP Polish**: ``±max(ulp_radius, 4)`` fp16
+     ULPs around both ``+d_sel`` and ``-d_sel``, scored on exact encode SSE.
+  5. **Closed-Form Vesica/DC Code Projection within ``sse_budget``**:
+     a bounded gradient project-and-accept step on the 3-bit codes using the
+     effective per-weight scale ``d*sc_j`` and the strict monotone
+     ``(SSE + spectral)`` decrease test.
+
+The final ``block_q3_K`` payload is emitted in the exact GGML layout
+(``hmask[32] + qs[64] + scales[12] + fp16 d``, 110 bytes).
 """
 
 from __future__ import annotations
@@ -40,6 +59,10 @@ QK3_K = 256
 Q3_K_BLOCK_BYTES = 110
 Q3_K_SUBBLOCK = 16
 Q3_K_N_SUBBLOCKS = QK3_K // Q3_K_SUBBLOCK
+Q3_SC_MIN = -32
+Q3_SC_MAX = 31
+Q3_CODE_MIN = -4
+Q3_CODE_MAX = 3
 Q3_N_CAND = 24
 Q3_N_STATES = 6
 
@@ -207,6 +230,11 @@ def _q3_sieve_and_viterbi(cand_errors: np.ndarray, cfg: HPCConfig) -> np.ndarray
     return np.where(g_e < cur_e * cfg.greedy_override_ratio, global_best, chosen_cand)
 
 
+# ---------------------------------------------------------------------------
+# Level A: HPC Q4_0 stack on 16-weight sub-blocks -> signed fp16 scales s_j
+# ---------------------------------------------------------------------------
+
+
 def _subblock_scales_hpc(
     X: np.ndarray, W: np.ndarray, cfg: HPCConfig
 ) -> np.ndarray:
@@ -318,91 +346,6 @@ def _subblock_scales_hpc(
     return best_d.astype(np.float32)
 
 
-def _fit_q3_k_blocks(
-    X: np.ndarray, W: np.ndarray, sub_scales: np.ndarray, cfg: HPCConfig
-) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
-    """Fit the GGML super-scale d and signed 6-bit scales for (B, 256) blocks."""
-    B = X.shape[0]
-    X = X.reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK).astype(np.float64)
-    W = W.reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK).astype(np.float64)
-    s = sub_scales.reshape(B, Q3_K_N_SUBBLOCKS).astype(np.float64)
-    rows = np.arange(B)
-
-    max_idx = np.argmax(np.abs(s), axis=1)
-    max_scale = s[rows, max_idx]
-    d0 = np.where(np.abs(max_scale) > 1.0e-30, -max_scale / 32.0, 0.0)
-
-    cand_d16 = to_fp16_bits((d0[:, None] * MULTS_D[None, :]).astype(np.float32))
-    cand_d = from_fp16_bits(cand_d16).astype(np.float64)
-    act_d = np.abs(cand_d) > 1.0e-15
-    safe_d = np.where(act_d, cand_d, 1.0)
-    sc = np.where(
-        act_d[:, :, None],
-        np.clip(nearest_int((s[:, None, :] / safe_d[:, :, None]).astype(np.float32)), -32, 31),
-        0,
-    )
-    scale = cand_d[:, :, None] * sc
-    act_s = np.abs(scale) > 1.0e-15
-    safe_s = np.where(act_s, scale, 1.0)
-    code = np.where(
-        act_s[:, :, :, None],
-        np.clip(nearest_int((X[:, None, :, :] / safe_s[:, :, :, None]).astype(np.float32)), -4, 3),
-        0,
-    )
-    e = X[:, None, :, :] - scale[:, :, :, None] * code
-    sse_c = np.sum(W[:, None, :, :] * e * e, axis=(2, 3))
-    best_cand = _q3_sieve_and_viterbi(sse_c, cfg)
-    best_d16 = cand_d16[rows, best_cand]
-
-    rad = max(int(cfg.ulp_radius), 4)
-    deltas = np.arange(-rad, rad + 1, dtype=np.int32)
-    pos_b = np.clip(
-        to_fp16_bits(np.abs(from_fp16_bits(best_d16))).astype(np.int32)[:, None] + deltas[None, :],
-        0,
-        0x7BFF,
-    ).astype(np.uint16)
-    neg_b = (pos_b | np.uint16(0x8000)).astype(np.uint16)
-    all_b = np.concatenate([pos_b, neg_b], axis=1)
-    all_d = from_fp16_bits(all_b).astype(np.float64)
-    act_u = np.abs(all_d) > 1.0e-15
-    safe_u = np.where(act_u, all_d, 1.0)
-    sc_u = np.where(
-        act_u[:, :, None],
-        np.clip(nearest_int((s[:, None, :] / safe_u[:, :, None]).astype(np.float32)), -32, 31),
-        0,
-    )
-    scale_u = all_d[:, :, None] * sc_u
-    act_su = np.abs(scale_u) > 1.0e-15
-    safe_su = np.where(act_su, scale_u, 1.0)
-    code_u = np.where(
-        act_su[:, :, :, None],
-        np.clip(nearest_int((X[:, None, :, :] / safe_su[:, :, :, None]).astype(np.float32)), -4, 3),
-        0,
-    )
-    e_u = X[:, None, :, :] - scale_u[:, :, :, None] * code_u
-    sse_u = np.sum(W[:, None, :, :] * e_u * e_u, axis=(2, 3))
-    pick_u = np.argmin(sse_u, axis=1)
-    best_d16 = all_b[rows, pick_u]
-    best_d = all_d[rows, pick_u]
-
-    act_b = np.abs(best_d) > 1.0e-15
-    safe_b = np.where(act_b, best_d, 1.0)
-    sc = np.where(
-        act_b[:, None],
-        np.clip(nearest_int((s / safe_b[:, None]).astype(np.float32)), -32, 31),
-        0,
-    )
-    scale = best_d[:, None] * sc
-    act_s = np.abs(scale) > 1.0e-15
-    safe_s = np.where(act_s, scale, 1.0)
-    code = np.where(
-        act_s[:, :, None],
-        np.clip(nearest_int((X / safe_s[:, :, None]).astype(np.float32)), -4, 3),
-        0,
-    )
-    return best_d16.astype(np.uint16), sc.astype(np.int16), code.astype(np.int16)
-
-
 def _subblock_scales_hpc_cuda(
     X_np: np.ndarray, W_np: np.ndarray, dev: "torch.device", cfg: HPCConfig
 ) -> np.ndarray:
@@ -522,6 +465,190 @@ def _subblock_scales_hpc_cuda(
         return best_d.cpu().numpy().astype(np.float32)
 
 
+# ---------------------------------------------------------------------------
+# Level B: GGML (d, sc, code) super-block fit — full Q4_0 decision stack
+# ---------------------------------------------------------------------------
+
+
+def _fit_q3_k_blocks(
+    X: np.ndarray, W: np.ndarray, sub_scales: np.ndarray, cfg: HPCConfig
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Full HPC V2 super-block fit for (B, 256) blocks on CPU."""
+    B = X.shape[0]
+    X3 = X.astype(np.float32).reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK)
+    W3 = W.astype(np.float32).reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK)
+    Xf = X3.reshape(B, QK3_K)
+    Wf = W3.reshape(B, QK3_K)
+    s = sub_scales.reshape(B, Q3_K_N_SUBBLOCKS).astype(np.float32)
+    rows = np.arange(B)
+    lam_dc = (float(cfg.dc_lambda) / QK3_K) if cfg.use_spectral else 0.0
+    lam_vw = (float(cfg.vw_lambda) / QK3_K) if cfg.use_spectral else 0.0
+    Wj = W3.sum(axis=2)  # (B, 16) per-sub-block importance
+
+    def _encode(d: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+        """JointLS quantize s with scale d: sc in [-32,31], codes in [-4,3]."""
+        act = np.abs(d) > 1.0e-15
+        safe_d = np.where(act, d, 1.0)
+        sc = np.where(
+            act[:, None],
+            np.clip(nearest_int((s / safe_d[:, None]).astype(np.float32)), Q3_SC_MIN, Q3_SC_MAX),
+            0,
+        )
+        scale = d[:, None] * sc.astype(np.float32)
+        act_s = np.abs(scale) > 1.0e-15
+        safe_s = np.where(act_s, scale, 1.0)
+        code = np.where(
+            act_s[:, :, None],
+            np.clip(nearest_int((X3 / safe_s[:, :, None]).astype(np.float32)), Q3_CODE_MIN, Q3_CODE_MAX),
+            0,
+        )
+        e = X3 - scale[:, :, None] * code.astype(np.float32)
+        sse = np.sum(W3 * e * e, axis=(1, 2), dtype=np.float64)
+        return sc, code, sse
+
+    # 1. Geometry-Preserving Analog + WLS Seeds for the super-scale d:
+    amax = np.max(np.abs(s), axis=1)
+    d_range = amax / 31.5
+    wsum = np.maximum(np.sum(Wj, axis=1, keepdims=True), 1.0e-30)
+    mu = np.sum(Wj * s, axis=1, keepdims=True) / wsum
+    var = np.sum(Wj * (s - mu) ** 2, axis=1, keepdims=True) / wsum
+    sig = np.sqrt(np.maximum(var, 1.0e-30))
+    s_comp = mu + sig * 1.85 * np.tanh((s - mu) / (1.85 * sig))
+    d_analog = (np.max(np.abs(s_comp), axis=1) / 31.5).astype(np.float32)
+
+    wls_d = d_range.copy()
+    best_seed_sse = np.full(B, np.inf, dtype=np.float64)
+    for d_start in (d_range, -d_range, d_analog, -d_analog):
+        d = d_start.astype(np.float64)
+        for _ in range(4):
+            act = np.abs(d) > 1.0e-15
+            sd = np.where(act, d, 1.0)[:, None]
+            sc = np.where(
+                act[:, None],
+                np.clip(nearest_int((s / sd).astype(np.float32)), Q3_SC_MIN, Q3_SC_MAX),
+                0,
+            )
+            c = sc.astype(np.float64)
+            num = np.sum(Wj * s * c, axis=1, dtype=np.float64)
+            den = np.sum(Wj * c * c, axis=1, dtype=np.float64)
+            valid = den > 1.0e-30
+            d = np.where(valid, num / np.where(valid, den, 1.0), 0.0)
+        df = from_fp16_bits(to_fp16_bits(d.astype(np.float32)))
+        _, _, sse = _encode(df)
+        take = sse < best_seed_sse
+        best_seed_sse = np.where(take, sse, best_seed_sse)
+        wls_d = np.where(take, df, wls_d)
+
+    # 2. 24-Candidate Super-Scale Lattice with Extended Spectral Objective:
+    cand_d16 = to_fp16_bits(wls_d[:, None] * MULTS_D[None, :])  # (B, 24)
+    cand_d = from_fp16_bits(cand_d16)
+    act_c = np.abs(cand_d) > 1.0e-15
+    safe_c = np.where(act_c, cand_d, 1.0)
+    sc_c = np.where(
+        act_c[:, :, None],
+        np.clip(nearest_int((s[:, None, :] / safe_c[:, :, None]).astype(np.float32)), Q3_SC_MIN, Q3_SC_MAX),
+        0,
+    )
+    scale_c = cand_d[:, :, None] * sc_c.astype(np.float32)  # (B, 24, 16)
+    act_s = np.abs(scale_c) > 1.0e-15
+    safe_s = np.where(act_s, scale_c, 1.0)
+    code_c = np.where(
+        act_s[:, :, :, None],
+        np.clip(nearest_int((X3[:, None] / safe_s[:, :, :, None]).astype(np.float32)), Q3_CODE_MIN, Q3_CODE_MAX),
+        0,
+    ).astype(np.int16)
+    e_f = (X3[:, None] - scale_c[:, :, :, None] * code_c.astype(np.float32)).reshape(B, Q3_N_CAND, QK3_K)
+    sse_c = np.sum(Wf[:, None, :] * e_f * e_f, axis=2)
+    if cfg.use_spectral and (lam_dc > 0.0 or lam_vw > 0.0):
+        half = QK3_K // 2
+        v_c = e_f[:, :, :half] + e_f[:, :, half:]
+        dc_c = np.sum(e_f, axis=2)
+        cand_errors = (sse_c + lam_vw * np.sum(v_c * v_c, axis=2) + lam_dc * dc_c * dc_c).astype(np.float32)
+    else:
+        cand_errors = sse_c.astype(np.float32)
+
+    # 3. 6-State Sieve + 6-State Viterbi DP over the super-scale lattice:
+    best_cand = _q3_sieve_and_viterbi(cand_errors, cfg)
+    d_sel = cand_d[rows, best_cand]
+
+    # 4. Asymmetric Sign-Flip fp16-ULP Polish (+/-8 ULPs around +d_sel and -d_sel):
+    rad = max(int(cfg.ulp_radius), 4)
+    deltas = np.arange(-rad, rad + 1, dtype=np.int32)
+    pos_b = np.clip(
+        to_fp16_bits(np.abs(d_sel)).astype(np.int32)[:, None] + deltas[None, :],
+        0,
+        0x7BFF,
+    ).astype(np.uint16)
+    neg_b = (pos_b | np.uint16(0x8000)).astype(np.uint16)
+    all_b = np.concatenate([pos_b, neg_b], axis=1)
+    all_d = from_fp16_bits(all_b)
+    act_u = np.abs(all_d) > 1.0e-15
+    safe_u = np.where(act_u, all_d, 1.0)
+    sc_u = np.where(
+        act_u[:, :, None],
+        np.clip(nearest_int((s[:, None, :] / safe_u[:, :, None]).astype(np.float32)), Q3_SC_MIN, Q3_SC_MAX),
+        0,
+    )
+    scale_u = all_d[:, :, None] * sc_u.astype(np.float32)
+    act_su = np.abs(scale_u) > 1.0e-15
+    safe_su = np.where(act_su, scale_u, 1.0)
+    code_u = np.where(
+        act_su[:, :, :, None],
+        np.clip(nearest_int((X3[:, None] / safe_su[:, :, :, None]).astype(np.float32)), Q3_CODE_MIN, Q3_CODE_MAX),
+        0,
+    )
+    e_u = X3[:, None] - scale_u[:, :, :, None] * code_u.astype(np.float32)
+    sse_u = np.sum(W3[:, None] * e_u * e_u, axis=(2, 3), dtype=np.float64)
+    pick_u = np.argmin(sse_u, axis=1)
+    best_d16 = all_b[rows, pick_u]
+    best_d = all_d[rows, pick_u]
+
+    # Final (sc, code) derived from the polished super-scale d:
+    act_b = np.abs(best_d) > 1.0e-15
+    safe_b = np.where(act_b, best_d, 1.0)
+    sc = np.where(
+        act_b[:, None],
+        np.clip(nearest_int((s / safe_b[:, None]).astype(np.float32)), Q3_SC_MIN, Q3_SC_MAX),
+        0,
+    )
+    scale = best_d[:, None] * sc.astype(np.float32)
+    act_s = np.abs(scale) > 1.0e-15
+    safe_s = np.where(act_s, scale, 1.0)
+    code = np.where(
+        act_s[:, :, None],
+        np.clip(nearest_int((X3 / safe_s[:, :, None]).astype(np.float32)), Q3_CODE_MIN, Q3_CODE_MAX),
+        0,
+    )
+
+    # 5. Closed-Form 1-D Quadratic Vesica + DC Code Projection within sse_budget:
+    if cfg.use_spectral and (lam_vw > 0.0 or lam_dc > 0.0):
+        d_w = np.repeat(best_d[:, None].astype(np.float64) * sc.astype(np.float64), Q3_K_SUBBLOCK, axis=1)
+        e = Xf.astype(np.float64) - d_w * code.reshape(B, QK3_K)
+        sse_live = np.sum(Wf * e * e, axis=1, dtype=np.float64)
+        sse_cap = sse_live * (1.0 + cfg.sse_budget)
+        half = QK3_K // 2
+        v = e[:, :half] + e[:, half:]
+        dc = np.sum(e, axis=1)
+        denom = (Wf + lam_vw + lam_dc) * (d_w * d_w)
+        act_bw = np.abs(d_w) > 1.0e-15
+        s_den = np.where(act_bw & (denom > 1.0e-30), denom, 1.0)
+        v_full = np.concatenate([v, v], axis=1)
+        grad = d_w * (Wf * e + lam_vw * v_full + lam_dc * dc[:, None])
+        dq = np.where(act_bw, nearest_int((grad / s_den).astype(np.float32)), 0)
+        code_f = code.reshape(B, QK3_K)
+        code_new = np.clip(code_f + dq, Q3_CODE_MIN, Q3_CODE_MAX)
+        e_new = Xf.astype(np.float64) - d_w * code_new
+        sse_new = np.sum(Wf * e_new * e_new, axis=1, dtype=np.float64)
+        v_new = e_new[:, :half] + e_new[:, half:]
+        dc_new = np.sum(e_new, axis=1)
+        spec_old = lam_vw * np.sum(v * v, axis=1) + lam_dc * dc * dc
+        spec_new = lam_vw * np.sum(v_new * v_new, axis=1) + lam_dc * dc_new * dc_new
+        accept = (sse_new <= sse_cap) & ((sse_new + spec_new) < (sse_live + spec_old) - 1.0e-15)
+        code = np.where(accept[:, None], code_new, code_f).reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK)
+
+    return best_d16.astype(np.uint16), sc.astype(np.int16), code.astype(np.int16)
+
+
 def _fit_q3_k_blocks_cuda(
     X_np: np.ndarray,
     W_np: np.ndarray,
@@ -529,74 +656,156 @@ def _fit_q3_k_blocks_cuda(
     dev: "torch.device",
     cfg: HPCConfig,
 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
-    """Fit GGML Q3_K super-scales and codes on CUDA."""
+    """Full HPC V2 super-block fit for (B, 256) blocks on CUDA."""
     B = X_np.shape[0]
+    lam_dc = (float(cfg.dc_lambda) / QK3_K) if cfg.use_spectral else 0.0
+    lam_vw = (float(cfg.vw_lambda) / QK3_K) if cfg.use_spectral else 0.0
+
     with torch.no_grad():
-        X = torch.from_numpy(X_np).to(dev, dtype=torch.float32).reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK)
-        W = torch.from_numpy(W_np).to(dev, dtype=torch.float32).reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK)
-        s = torch.from_numpy(sub_scales.reshape(B, Q3_K_N_SUBBLOCKS)).to(dev, dtype=torch.float32)
+        X3 = torch.from_numpy(X_np).to(dev, dtype=torch.float32).reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK)
+        W3 = torch.from_numpy(W_np).to(dev, dtype=torch.float32).reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK)
+        Xf = X3.reshape(B, QK3_K)
+        Wf = W3.reshape(B, QK3_K)
+        s = torch.from_numpy(np.ascontiguousarray(sub_scales.reshape(B, Q3_K_N_SUBBLOCKS))).to(dev, dtype=torch.float32)
         rows = torch.arange(B, device=dev)
+        Wj = W3.sum(dim=2)  # (B, 16) per-sub-block importance
+
+        def _encode(d: "torch.Tensor"):
+            """JointLS quantize s with scale d: sc in [-32,31], codes in [-4,3]."""
+            act = torch.abs(d) > 1.0e-15
+            safe_d = torch.where(act, d, torch.ones_like(d))
+            q_sc = _nearest_int_t(s / safe_d[:, None]).clamp(Q3_SC_MIN, Q3_SC_MAX)
+            sc = torch.where(act[:, None], q_sc, torch.zeros_like(q_sc))
+            scale = d[:, None] * sc.to(torch.float32)
+            act_s = torch.abs(scale) > 1.0e-15
+            safe_s = torch.where(act_s, scale, torch.ones_like(scale))
+            q_cd = _nearest_int_t(X3 / safe_s[:, :, None]).clamp(Q3_CODE_MIN, Q3_CODE_MAX)
+            code = torch.where(act_s[:, :, None], q_cd, torch.zeros_like(q_cd))
+            e = X3 - scale[:, :, None] * code.to(torch.float32)
+            sse = torch.sum(W3 * e * e, dim=(1, 2))
+            return sc, code, sse
+
+        # 1. Geometry-Preserving Analog + WLS Seeds for the super-scale d:
+        amax = torch.max(torch.abs(s), dim=1).values
+        d_range = amax / 31.5
+        wsum = torch.clamp(torch.sum(Wj, dim=1, keepdim=True), min=1.0e-30)
+        mu = torch.sum(Wj * s, dim=1, keepdim=True) / wsum
+        var = torch.sum(Wj * (s - mu) ** 2, dim=1, keepdim=True) / wsum
+        sig = torch.sqrt(torch.clamp(var, min=1.0e-30))
+        s_comp = mu + sig * 1.85 * torch.tanh((s - mu) / (1.85 * sig))
+        d_analog = torch.max(torch.abs(s_comp), dim=1).values / 31.5
+
+        wls_d = d_range.clone()
+        best_seed_sse = torch.full((B,), float("inf"), dtype=torch.float32, device=dev)
+        for d_start in (d_range, -d_range, d_analog, -d_analog):
+            d = d_start.clone()
+            for _ in range(4):
+                act = torch.abs(d) > 1.0e-15
+                sd = torch.where(act, d, torch.ones_like(d))[:, None]
+                qi = _nearest_int_t(s / sd).clamp(Q3_SC_MIN, Q3_SC_MAX)
+                sc = torch.where(act[:, None], qi, torch.zeros_like(qi))
+                c = sc.to(torch.float32)
+                num = torch.sum(Wj * s * c, dim=1)
+                den = torch.sum(Wj * c * c, dim=1)
+                valid = den > 1.0e-30
+                d = torch.where(valid, num / torch.where(valid, den, torch.ones_like(den)), torch.zeros_like(den))
+            df = d.to(torch.float16).to(torch.float32)
+            _, _, sse = _encode(df)
+            take = sse < best_seed_sse
+            best_seed_sse = torch.where(take, sse, best_seed_sse)
+            wls_d = torch.where(take, df, wls_d)
 
-        max_idx = torch.argmax(torch.abs(s), dim=1)
-        max_scale = s[rows, max_idx]
-        d0 = torch.where(torch.abs(max_scale) > 1.0e-30, -max_scale / 32.0, torch.zeros_like(max_scale))
+        # 2. 24-Candidate Super-Scale Lattice with Extended Spectral Objective:
         mults = torch.tensor(MULTS_D, dtype=torch.float32, device=dev)
-        cand_d = (d0[:, None] * mults[None, :]).to(torch.float16).to(torch.float32)
-        cand_bits = cand_d.to(torch.float16).view(torch.int16)
-
-        act_d = torch.abs(cand_d) > 1.0e-15
-        safe_d = torch.where(act_d, cand_d, torch.ones_like(cand_d))
-        qi = _nearest_int_t(s[:, None, :] / safe_d[:, :, None])
-        sc = torch.where(act_d[:, :, None], torch.clamp(qi, -32, 31), torch.zeros_like(qi))
-        scale = cand_d[:, :, None] * sc.to(torch.float32)
-        act_s = torch.abs(scale) > 1.0e-15
-        safe_s = torch.where(act_s, scale, torch.ones_like(scale))
-        cq = _nearest_int_t(X[:, None, :, :] / safe_s[:, :, :, None])
-        code = torch.where(act_s[:, :, :, None], torch.clamp(cq, -4, 3), torch.zeros_like(cq))
-        e = X[:, None, :, :] - scale[:, :, :, None] * code.to(torch.float32)
-        sse_c = torch.sum(W[:, None, :, :] * e * e, dim=(2, 3))
+        cand_d = (wls_d[:, None] * mults[None, :]).to(torch.float16).to(torch.float32)
+        act_c = torch.abs(cand_d) > 1.0e-15
+        safe_c = torch.where(act_c, cand_d, torch.ones_like(cand_d))
+        q_sc_c = _nearest_int_t(s[:, None, :] / safe_c[:, :, None]).clamp(Q3_SC_MIN, Q3_SC_MAX)
+        sc_c = torch.where(act_c[:, :, None], q_sc_c, torch.zeros_like(q_sc_c))
+        scale_c = cand_d[:, :, None] * sc_c.to(torch.float32)
+        act_s = torch.abs(scale_c) > 1.0e-15
+        safe_s = torch.where(act_s, scale_c, torch.ones_like(scale_c))
+        q_cd_c = _nearest_int_t(X3[:, None] / safe_s[:, :, :, None]).clamp(Q3_CODE_MIN, Q3_CODE_MAX)
+        code_c = torch.where(act_s[:, :, :, None], q_cd_c, torch.zeros_like(q_cd_c))
+        e_f = (X3[:, None] - scale_c[:, :, :, None] * code_c.to(torch.float32)).reshape(B, Q3_N_CAND, QK3_K)
+        sse_c = torch.sum(Wf[:, None, :] * e_f * e_f, dim=2)
+        if cfg.use_spectral and (lam_dc > 0.0 or lam_vw > 0.0):
+            half = QK3_K // 2
+            v_c = e_f[:, :, :half] + e_f[:, :, half:]
+            dc_c = torch.sum(e_f, dim=2)
+            cand_err_t = sse_c + lam_vw * torch.sum(v_c * v_c, dim=2) + lam_dc * dc_c * dc_c
+        else:
+            cand_err_t = sse_c
 
-        best_cand_np = _q3_sieve_and_viterbi(sse_c.cpu().numpy(), cfg)
+        # 3. 6-State Sieve + 6-State Viterbi DP over the super-scale lattice:
+        best_cand_np = _q3_sieve_and_viterbi(cand_err_t.cpu().numpy(), cfg)
         best_cand_t = torch.from_numpy(best_cand_np).to(dev)
-        best_bits = cand_bits[rows, best_cand_t]
+        d_sel = cand_d[rows, best_cand_t]
 
+        # 4. Asymmetric Sign-Flip fp16-ULP Polish (+/-8 ULPs around +d_sel and -d_sel):
         rad = max(int(cfg.ulp_radius), 4)
         deltas = torch.arange(-rad, rad + 1, dtype=torch.int32, device=dev)
-        abs_bits = best_bits.to(torch.int32) & 0x7FFF
+        abs_bits = torch.abs(d_sel).to(torch.float16).view(torch.int16).to(torch.int32) & 0x7FFF
         pos_b = torch.clamp(abs_bits[:, None] + deltas[None, :], 0, 0x7BFF)
         neg_b = pos_b | 0x8000
         all_b = torch.cat([pos_b, neg_b], dim=1)
         all_d = all_b.to(torch.int16).view(torch.float16).to(torch.float32)
         act_u = torch.abs(all_d) > 1.0e-15
         safe_u = torch.where(act_u, all_d, torch.ones_like(all_d))
-        q_u = _nearest_int_t(s[:, None, :] / safe_u[:, :, None])
-        sc_u = torch.where(act_u[:, :, None], torch.clamp(q_u, -32, 31), torch.zeros_like(q_u))
+        q_sc_u = _nearest_int_t(s[:, None, :] / safe_u[:, :, None]).clamp(Q3_SC_MIN, Q3_SC_MAX)
+        sc_u = torch.where(act_u[:, :, None], q_sc_u, torch.zeros_like(q_sc_u))
         scale_u = all_d[:, :, None] * sc_u.to(torch.float32)
         act_su = torch.abs(scale_u) > 1.0e-15
         safe_su = torch.where(act_su, scale_u, torch.ones_like(scale_u))
-        c_u = _nearest_int_t(X[:, None, :, :] / safe_su[:, :, :, None])
-        code_u = torch.where(act_su[:, :, :, None], torch.clamp(c_u, -4, 3), torch.zeros_like(c_u))
-        e_u = X[:, None, :, :] - scale_u[:, :, :, None] * code_u.to(torch.float32)
-        sse_u = torch.sum(W[:, None, :, :] * e_u * e_u, dim=(2, 3))
+        q_cd_u = _nearest_int_t(X3[:, None] / safe_su[:, :, :, None]).clamp(Q3_CODE_MIN, Q3_CODE_MAX)
+        code_u = torch.where(act_su[:, :, :, None], q_cd_u, torch.zeros_like(q_cd_u))
+        e_u = X3[:, None] - scale_u[:, :, :, None] * code_u.to(torch.float32)
+        sse_u = torch.sum(W3[:, None] * e_u * e_u, dim=(2, 3))
         pick_u = torch.argmin(sse_u, dim=1)
         best_d16 = all_b[rows, pick_u]
         best_d = all_d[rows, pick_u]
 
+        # Final (sc, code) derived from the polished super-scale d:
         act_b = torch.abs(best_d) > 1.0e-15
         safe_b = torch.where(act_b, best_d, torch.ones_like(best_d))
-        sc = torch.where(
-            act_b[:, None],
-            torch.clamp(_nearest_int_t(s / safe_b[:, None]), -32, 31),
-            torch.zeros_like(s, dtype=torch.int32),
-        )
+        q_bs = _nearest_int_t(s / safe_b[:, None]).clamp(Q3_SC_MIN, Q3_SC_MAX)
+        sc = torch.where(act_b[:, None], q_bs, torch.zeros_like(q_bs))
         scale = best_d[:, None] * sc.to(torch.float32)
         act_s = torch.abs(scale) > 1.0e-15
         safe_s = torch.where(act_s, scale, torch.ones_like(scale))
-        code = torch.where(
-            act_s[:, :, None],
-            torch.clamp(_nearest_int_t(X / safe_s[:, :, None]), -4, 3),
-            torch.zeros_like(sc, dtype=torch.int32)[:, :, None].expand(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK),
-        )
+        q_bd = _nearest_int_t(X3 / safe_s[:, :, None]).clamp(Q3_CODE_MIN, Q3_CODE_MAX)
+        code = torch.where(act_s[:, :, None], q_bd, torch.zeros_like(q_bd))
+        code_f = code.reshape(B, QK3_K)
+
+        # 5. Closed-Form Vesica/DC Code Projection within sse_budget:
+        if cfg.use_spectral and (lam_vw > 0.0 or lam_dc > 0.0):
+            d_w = (best_d[:, None] * sc.to(torch.float32)).repeat_interleave(Q3_K_SUBBLOCK, dim=1)
+            e_f64 = Xf - d_w * code_f.to(torch.float32)
+            sse_live = torch.sum(Wf * e_f64 * e_f64, dim=1)
+            sse_cap = sse_live * (1.0 + float(cfg.sse_budget))
+            half = QK3_K // 2
+            v = e_f64[:, :half] + e_f64[:, half:]
+            dc = torch.sum(e_f64, dim=1)
+            denom = (Wf + lam_vw + lam_dc) * (d_w * d_w)
+            act_bw = torch.abs(d_w) > 1.0e-15
+            s_den = torch.where(act_bw & (denom > 1.0e-30), denom, torch.ones_like(denom))
+            v_full = torch.cat([v, v], dim=1)
+            grad = d_w * (Wf * e_f64 + lam_vw * v_full + lam_dc * dc[:, None])
+            dq = torch.where(
+                act_bw,
+                _nearest_int_t(grad / s_den),
+                torch.zeros_like(code_f),
+            )
+            code_new = torch.clamp(code_f + dq, Q3_CODE_MIN, Q3_CODE_MAX)
+            e_new = Xf - d_w * code_new.to(torch.float32)
+            sse_new = torch.sum(Wf * e_new * e_new, dim=1)
+            v_new = e_new[:, :half] + e_new[:, half:]
+            dc_new = torch.sum(e_new, dim=1)
+            spec_old = lam_vw * torch.sum(v * v, dim=1) + lam_dc * dc * dc
+            spec_new = lam_vw * torch.sum(v_new * v_new, dim=1) + lam_dc * dc_new * dc_new
+            accept = (sse_new <= sse_cap) & ((sse_new + spec_new) < (sse_live + spec_old) - 1.0e-12)
+            code = torch.where(accept[:, None], code_new, code_f).reshape(B, Q3_K_N_SUBBLOCKS, Q3_K_SUBBLOCK)
+
         return (
             best_d16.cpu().numpy().astype(np.uint16),
             sc.cpu().numpy().astype(np.int16),
@@ -685,7 +894,7 @@ def quantize_tensor_q3_k(
     else:
         W_raw = None
 
-    use_cuda = cfg.device != "cpu" and torch is not None and cuda_is_available() and (
+    use_cuda = cfg.device != "cpu" and cuda_is_available() and (
         cfg.device.startswith("cuda") or (cfg.device == "auto" and n_blocks >= 64)
     )
     dev = torch.device("cuda" if cfg.device == "auto" else cfg.device) if use_cuda else None