diff --git a/NEWS.md b/NEWS.md
index 72b2fbab..aaff3adf 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -50,6 +50,27 @@
   report CL2 leverage diagnostics: the returned covariance carries
   `max_leverage` and `n_capped_clusters` attributes, and a warning is emitted
   when one or more clusters hit the leverage cap.
+* CL2 now uses the exact Bell--McCaffrey leverage block by default when the
+  finite-sample correction is relevant and the dense-block calculation is
+  affordable (at most 100 clusters and cost proxy `G * max(D_g) * p^2 <=
+  5e9`, a measured sub-second marginal cost). At larger/prohibitive problems it uses the historical shortcut.
+  `cl2_adjustment = "exact"` or `"shortcut"` on `pffr()` and `coef.pffr()`
+  explicitly selects either variant. Exact-CL2 diagnostics include the number
+  of eigenvalue-floored blocks (`n_adjusted`).
+* The CL2 sandwich now checks the penalized hat matrix against its own bounds
+  (`0 <= h_ii <= 1`, `0 <= eigen(H_gg) <= 1`, and positive semi-definiteness of
+  the exact Bell--McCaffrey block) and **warns** when they are violated by more
+  than round-off. Such a violation means the model-based bread and the weighted
+  design have become numerically inconsistent -- a barely converged or
+  extremely ill-conditioned fit -- so the cluster-robust covariance is
+  meaningless however plausible it looks. Previously this produced silently
+  exploded standard errors (interval widths up to 1e133 were observed on
+  degenerate Poisson fits, with nothing to distinguish them from a legitimately
+  wide interval). The covariance is still returned, now carrying
+  `max_obs_leverage` and `hat_invariant_violation` attributes. The CR1
+  (`sandwich = "cluster"`) covariance is built from the same bread and is
+  equally affected, so switching sandwich type is not a remedy: inspect and
+  refit the model.
 * AR(1) support improvements: `pffr()` now automatically switches to
   `algorithm = "bam"` and `method = "fREML"` when `rho` is supplied, and
   sets `discrete = TRUE` for non-Gaussian families.
diff --git a/R/pffr-core.R b/R/pffr-core.R
index 61c8567a..354ab921 100644
--- a/R/pffr-core.R
+++ b/R/pffr-core.R
@@ -675,7 +675,8 @@ pffr_build_metadata <- function(
   ydata,
   sandwich,
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 ) {
   list(
     call = call,
@@ -698,6 +699,7 @@ pffr_build_metadata <- function(
     sandwich = sandwich,
     dof_correction = dof_correction,
     edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment,
     # Covariance storage contract version: format 2 keeps $Vp/$Vc/$Ve
     # model-based ALWAYS; the robust covariance lives in $pffr$Vsandwich.
     cov_format = PFFR_COV_STORAGE_FORMAT,
@@ -1616,6 +1618,10 @@ gam_sandwich_cluster <- function(
 #'   If `FALSE` (default), use Bayesian sandwich (`B2 = Vp - Ve`).
 #' @param tol Eigenvalue floor for numerical stability.
 #' @param leverage_cap Cap for cluster leverage eigenvalues (< 1).
+#' @param cl2_adjustment CL2 leverage adjustment: `"auto"` (default) selects
+#'   the exact Bell--McCaffrey block where it is relevant and affordable,
+#'   `"exact"` forces that block, and `"shortcut"` uses the historical
+#'   per-cluster shortcut `(I - H_gg)^(-1/2)`.
 #' @param b2 Internal ablation switch (default `TRUE` = current behavior). When
 #'   `FALSE`, drop the additive Bayesian smoothing-bias term \eqn{B_2 = V_p -
 #'   V_e} (X5).
@@ -1623,8 +1629,13 @@ gam_sandwich_cluster <- function(
 #'   behavior). When `TRUE`, center the leverage-adjusted per-cluster
 #'   contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat
 #'   (X6).
-#' @returns A p x p covariance matrix with attributes `n_capped_clusters` and
-#'   `max_leverage` for the CL2 leverage diagnostic.
+#' @returns A p x p covariance matrix with leverage diagnostics. Exact CL2
+#'   returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`),
+#'   `min_block_eig`, `max_block_kappa` and `min_block_eig_rel`; the shortcut
+#'   returns the legacy `n_capped_clusters` and `max_leverage` attributes. Both
+#'   return `max_obs_leverage` and `hat_invariant_violation` (`NULL` when the
+#'   penalized hat respects its bounds, otherwise a description of the
+#'   violation; see [pffr_hat_invariant_violation()]).
 #' @keywords internal
 gam_sandwich_cluster_cl2 <- function(
   b,
@@ -1632,6 +1643,7 @@ gam_sandwich_cluster_cl2 <- function(
   freq = FALSE,
   tol = 1e-8,
   leverage_cap = 0.999,
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   b2 = TRUE,
   center_scores = FALSE
 ) {
@@ -1673,41 +1685,98 @@ gam_sandwich_cluster_cl2 <- function(
   G <- n_clusters_checked(cluster_id_work)
   groups <- unique(cluster_id_work)
 
-  Vp <- b$Vp
-  B2 <- if (freq) 0 else b$Vp - b$Ve
+  Vp_raw <- b$Vp
+  Vp <- 0.5 * (Vp_raw + t(Vp_raw))
   p <- ncol(Xw)
+  cl2_adjustment <- resolve_cl2_adjustment(
+    cl2_adjustment,
+    G = G,
+    maxDg = max(table(cluster_id_work)),
+    p = p
+  )
+  B2 <- if (freq) {
+    0
+  } else if (cl2_adjustment == "exact") {
+    Vp - b$Ve
+  } else {
+    Vp_raw - b$Ve
+  }
+  Cmat <- if (cl2_adjustment == "exact") crossprod(Xw) else NULL
   meat <- matrix(0, nrow = p, ncol = p)
   Usum <- numeric(p)
   n_capped_clusters <- 0L
   max_leverage <- NA_real_
+  n_adjusted <- 0L
+  min_block_eig <- Inf
+  max_block_kappa <- 0
+  # Hat-invariant monitors (see pffr_hat_invariant_violation()). Both are free:
+  # h_ii is a by-product of T_g, and the block eigenvalues are computed anyway.
+  max_obs_leverage <- NA_real_
+  min_block_eig_rel <- 0
 
   for (g in groups) {
     idx <- which(cluster_id_work == g)
     Xwg <- Xw[idx, , drop = FALSE]
     zg <- z[idx]
 
-    Hgg <- Xwg %*% Vp %*% t(Xwg)
+    Tg <- Xwg %*% Vp
+    Hgg <- Tg %*% t(Xwg)
     Hgg <- 0.5 * (Hgg + t(Hgg))
 
-    ee_H <- eigen(Hgg, symmetric = TRUE)
-    cluster_max_leverage <- max(ee_H$values, na.rm = TRUE)
-    if (is.finite(cluster_max_leverage)) {
-      max_leverage <- if (is.na(max_leverage)) {
-        cluster_max_leverage
-      } else {
-        max(max_leverage, cluster_max_leverage)
-      }
-    }
-    if (any(ee_H$values > leverage_cap, na.rm = TRUE)) {
-      n_capped_clusters <- n_capped_clusters + 1L
-      ee_H$values <- pmin(ee_H$values, leverage_cap)
-      Hgg <- ee_H$vectors %*%
-        diag(ee_H$values, nrow = length(ee_H$values)) %*%
-        t(ee_H$vectors)
+    # Per-observation leverage h_ii = diag(H)_ii, mathematically in [0, 1] for
+    # the penalized hat. Tracked in both branches; a value above 1 means the
+    # bread and the weighted design are numerically inconsistent.
+    h_ii <- rowSums(Tg * Xwg)
+    if (any(is.finite(h_ii))) {
+      max_obs_leverage <- max(max_obs_leverage, h_ii, na.rm = TRUE)
     }
 
-    Mg <- diag(length(idx)) - Hgg
-    Ag <- sym_inv_sqrt(Mg, tol = tol)
+    if (cl2_adjustment == "shortcut") {
+      ee_H <- eigen(Hgg, symmetric = TRUE)
+      cluster_max_leverage <- max(ee_H$values, na.rm = TRUE)
+      if (is.finite(cluster_max_leverage)) {
+        max_leverage <- if (is.na(max_leverage)) {
+          cluster_max_leverage
+        } else {
+          max(max_leverage, cluster_max_leverage)
+        }
+      }
+      if (any(ee_H$values > leverage_cap, na.rm = TRUE)) {
+        n_capped_clusters <- n_capped_clusters + 1L
+        ee_H$values <- pmin(ee_H$values, leverage_cap)
+        Hgg <- ee_H$vectors %*%
+          diag(ee_H$values, nrow = length(ee_H$values)) %*%
+          t(ee_H$vectors)
+      }
+      Ag <- sym_inv_sqrt(diag(length(idx)) - Hgg, tol = tol)
+    } else {
+      # The exact Bell--McCaffrey block is B_g = I - 2 H_gg +
+      # (H_t^2)_gg, with (H_t^2)_gg = (Xw_g Vp) (Xw' Xw)
+      # (Xw_g Vp)'. This is a squared operator, so the gain-matched analogue
+      # of capping H_gg at `leverage_cap` floors B_g at (1 - cap)^2.
+      Bg <- diag(length(idx)) - 2 * Hgg + Tg %*% Cmat %*% t(Tg)
+      ee_B <- eigen(0.5 * (Bg + t(Bg)), symmetric = TRUE)
+      block_min <- min(ee_B$values)
+      min_block_eig <- min(min_block_eig, block_min)
+      # B_g = ((I - H)^2)_gg is a principal block of a squared symmetric
+      # matrix and is therefore positive semi-definite in exact arithmetic.
+      # A materially negative eigenvalue, measured relative to the block's own
+      # scale, is an invariant violation rather than round-off.
+      min_block_eig_rel <- min(
+        min_block_eig_rel,
+        block_min / max(abs(max(ee_B$values)), 1e-300)
+      )
+      max_block_kappa <- max(
+        max_block_kappa,
+        abs(max(ee_B$values)) / max(abs(block_min), 1e-300)
+      )
+      floor_val <- (1 - leverage_cap)^2
+      if (block_min < floor_val) n_adjusted <- n_adjusted + 1L
+      ee_B$values <- pmax(ee_B$values, floor_val)
+      Ag <- ee_B$vectors %*%
+        diag(1 / sqrt(ee_B$values), nrow = length(ee_B$values)) %*%
+        t(ee_B$vectors)
+    }
     Ug <- crossprod(Xwg, Ag %*% zg)
     meat <- meat + Ug %*% t(Ug)
     Usum <- Usum + as.vector(Ug)
@@ -1723,9 +1792,38 @@ gam_sandwich_cluster_cl2 <- function(
   V <- hc1 * Vp %*% meat %*% Vp
   if (isTRUE(b2)) V <- V + B2
   V <- 0.5 * (V + t(V))
+  attr(V, "cl2_adjustment") <- cl2_adjustment
   attr(V, "n_capped_clusters") <- n_capped_clusters
   attr(V, "max_leverage") <- max_leverage
-  if (n_capped_clusters > 0) {
+  attr(V, "n_adjusted") <- n_adjusted
+  attr(V, "min_block_eig") <- if (is.finite(min_block_eig)) min_block_eig else
+    NA_real_
+  attr(V, "max_block_kappa") <- max_block_kappa
+  attr(V, "max_obs_leverage") <- max_obs_leverage
+  attr(V, "min_block_eig_rel") <- min_block_eig_rel
+
+  hat_violation <- pffr_hat_invariant_violation(
+    max_obs_leverage = max_obs_leverage,
+    max_leverage = if (cl2_adjustment == "shortcut") max_leverage else NA_real_,
+    min_block_eig_rel = if (cl2_adjustment == "exact") min_block_eig_rel else 0
+  )
+  attr(V, "hat_invariant_violation") <- hat_violation
+  if (!is.null(hat_violation)) {
+    warning(
+      "Cluster-robust covariance is NOT trustworthy for this fit: ",
+      hat_violation,
+      " The penalized hat matrix satisfies 0 <= h_ii <= 1 and ",
+      "0 <= eigen(H_gg) <= 1 exactly, so this indicates that the model-based ",
+      "bread and the weighted design have become numerically inconsistent -- ",
+      "typically an ill-conditioned or barely converged fit (extreme fitted ",
+      "values, a near-singular penalized Hessian, or a basis far too rich for ",
+      "the data). Both the CL2 and the CR1 cluster-robust covariances inherit ",
+      "the problem, so switching sandwich type does not repair it: inspect and ",
+      "refit the model instead.",
+      call. = FALSE
+    )
+  }
+  if (cl2_adjustment == "shortcut" && n_capped_clusters > 0) {
     max_leverage_label <- if (is.finite(max_leverage)) {
       sprintf("%.3f", max_leverage)
     } else {
@@ -1750,6 +1848,121 @@ gam_sandwich_cluster_cl2 <- function(
   V
 }
 
+#' Numerical-sanity check on the per-cluster hat blocks
+#'
+#' The penalized hat \eqn{H = X_w V_p X_w'} obeys \eqn{0 \le h_{ii} \le 1} and,
+#' since any principal block of a symmetric matrix has its eigenvalues inside
+#' the parent's range, \eqn{0 \le \mathrm{eigen}(H_{gg}) \le 1}. The exact
+#' Bell--McCaffrey block \eqn{B_g = ((I - H)^2)_{gg}} is positive semi-definite
+#' for the same reason. When a fit is numerically degenerate --- an
+#' ill-conditioned or barely converged fit, e.g. a Poisson fit whose fitted
+#' means span many orders of magnitude --- the bread \eqn{V_p} stops being the
+#' inverse of the same weighted cross-product, and these invariants break by
+#' orders of magnitude rather than by round-off. This was the mechanism behind
+#' the numerically exploded interval widths observed on hard Poisson
+#' replicates in the locked benchmark (study LB, claim P-LB5).
+#'
+#' Detection is deliberately separate from the leverage cap. Capping
+#' \eqn{H_{gg}} at `leverage_cap` silently converts an impossible leverage into
+#' the *largest legitimate* one, so the shortcut adjustment responds to garbage
+#' with its maximal variance inflation; the exact block, being a squared
+#' operator, instead deflates. Both are bounded, neither is meaningful, and
+#' only a diagnostic can tell the user which situation they are in. We
+#' therefore warn and keep going rather than substituting a differently wrong
+#' number: the CR1 covariance is built from the same bread and is no more
+#' trustworthy here.
+#'
+#' @param max_obs_leverage Largest per-observation leverage \eqn{h_{ii}} seen,
+#'   or `NA`.
+#' @param max_leverage Largest eigenvalue of any \eqn{H_{gg}} (shortcut path
+#'   only; `NA` otherwise).
+#' @param min_block_eig_rel Smallest eigenvalue of any \eqn{B_g} relative to
+#'   that block's largest eigenvalue (exact path only; `0` otherwise).
+#' @param tol Relative slack allowed before an invariant counts as violated.
+#' @returns `NULL` when every invariant holds, otherwise a one-sentence
+#'   character description of the violation.
+#' @keywords internal
+pffr_hat_invariant_violation <- function(
+  max_obs_leverage = NA_real_,
+  max_leverage = NA_real_,
+  min_block_eig_rel = 0,
+  tol = 1e-6
+) {
+  msgs <- character(0)
+  if (isTRUE(is.finite(max_obs_leverage) && max_obs_leverage > 1 + tol)) {
+    msgs <- c(
+      msgs,
+      sprintf(
+        "the largest per-observation leverage is %.3g, above the bound 1;",
+        max_obs_leverage
+      )
+    )
+  }
+  if (isTRUE(is.finite(max_leverage) && max_leverage > 1 + tol)) {
+    msgs <- c(
+      msgs,
+      sprintf(
+        "the largest per-cluster hat eigenvalue is %.3g, above the bound 1;",
+        max_leverage
+      )
+    )
+  }
+  if (isTRUE(is.finite(min_block_eig_rel) && min_block_eig_rel < -tol)) {
+    msgs <- c(
+      msgs,
+      sprintf(
+        paste0(
+          "the exact Bell-McCaffrey block has a relative eigenvalue of %.3g, ",
+          "but it is positive semi-definite by construction;"
+        ),
+        min_block_eig_rel
+      )
+    )
+  }
+  if (length(msgs) == 0) return(NULL)
+  paste(msgs, collapse = " ")
+}
+
+#' Resolve the CL2 leverage adjustment
+#'
+#' Exact CL2 is the default when the finite-sample correction can matter
+#' (`G <= 100`; Study EX measured the exact-over-shortcut coverage gain
+#' decaying from ~+1pp at G = 20 to under +0.15pp at G = 100) and the dense
+#' block multiplication is affordable. The cost proxy `G * maxDg * p^2`
+#' covers the dominant `(Xw_g Vp) (Xw'Xw)` multiplication. The 5e9 bound is
+#' calibrated by measurement, not guessed: both adjustments pay the same
+#' per-cluster eigendecomposition, so the exact block's MARGINAL cost is only
+#' these multiplications -- timed at ~0.03 s extra at proxy 5e8 and ~0.7 s
+#' total at 5e9 (single-thread BLAS), i.e. the bound caps the correctness
+#' upgrade at well under a second and in particular never denies it in the
+#' saturated small-G regime where its benefit is largest. Large-G CL2 uses
+#' the historical shortcut because the correction tends to zero there.
+#'
+#' @param cl2_adjustment One of `"auto"`, `"exact"`, or `"shortcut"`.
+#' @param G Number of clusters.
+#' @param maxDg Largest cluster size.
+#' @param p Coefficient-space dimension.
+#' @returns `"exact"` or `"shortcut"`.
+#' @keywords internal
+resolve_cl2_adjustment <- function(
+  cl2_adjustment = c("auto", "exact", "shortcut"),
+  G,
+  maxDg,
+  p
+) {
+  cl2_adjustment <- match.arg(cl2_adjustment)
+  if (cl2_adjustment != "auto") return(cl2_adjustment)
+
+  cost <- G * maxDg * p^2
+  if (
+    is.finite(G) && is.finite(maxDg) && is.finite(p) && G <= 100 && cost <= 5e9
+  ) {
+    "exact"
+  } else {
+    "shortcut"
+  }
+}
+
 #' Resolve a pffr fit's covariance matrices to a canonical representation
 #'
 #' The current storage contract (format 2) keeps `Vp`/`Vc`/`Ve` model-based
@@ -1899,6 +2112,8 @@ restore_model_cov <- function(object) {
 #' @param center_scores Internal ablation switch (default `FALSE`); center the
 #'   per-cluster score sums before the meat when `TRUE` (X6). Applies to
 #'   `"cluster"`/`"cl2"` only.
+#' @param cl2_adjustment CL2 leverage adjustment (`"auto"`, `"exact"`, or
+#'   `"shortcut"`), used only for `type = "cl2"`.
 #' @returns A covariance matrix (with CL2 leverage attributes for `type =
 #'   "cl2"`).
 #' @keywords internal
@@ -1910,7 +2125,8 @@ pffr_compute_sandwich <- function(
   dof_correction = "none",
   edf_type = "trace",
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = "auto"
 ) {
   switch(
     type,
@@ -1928,7 +2144,8 @@ pffr_compute_sandwich <- function(
       cluster_id,
       freq = freq,
       b2 = b2,
-      center_scores = center_scores
+      center_scores = center_scores,
+      cl2_adjustment = cl2_adjustment
     ),
     hc = mgcv::vcov.gam(b, sandwich = TRUE, freq = freq),
     none = if (freq) b$Ve else (b$Vc %||% b$Vp),
@@ -1960,6 +2177,8 @@ pffr_compute_sandwich <- function(
 #' @param center_scores Internal ablation switch (default `FALSE` = current
 #'   behavior). When `TRUE`, center the per-cluster score sums before forming
 #'   the cluster/CL2 meat (X6). Same cache semantics as `b2`.
+#' @param cl2_adjustment CL2 leverage adjustment. `NULL` (default) inherits
+#'   the fitted choice; otherwise one of `"auto"`, `"exact"`, or `"shortcut"`.
 #' @returns A covariance matrix.
 #' @keywords internal
 pffr_vcov <- function(
@@ -1970,7 +2189,8 @@ pffr_vcov <- function(
   dof_correction = NULL,
   edf_type = NULL,
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = NULL
 ) {
   # Ablation variants (X5/X6) bypass the fit-time and recompute caches entirely,
   # so they never overwrite or shadow the standard cached matrices.
@@ -1990,6 +2210,13 @@ pffr_vcov <- function(
 
   dof_correction <- dof_correction %||% (object$pffr$dof_correction %||% "none")
   edf_type <- edf_type %||% (object$pffr$edf_type %||% "trace")
+  cl2_adjustment_explicit <- !is.null(cl2_adjustment)
+  cl2_adjustment <- cl2_adjustment %||%
+    (object$pffr$cl2_adjustment %||% "auto")
+  cl2_adjustment <- match.arg(
+    cl2_adjustment,
+    c("auto", "exact", "shortcut")
+  )
 
   # Serve the cached fit-time robust matrix when the request matches it exactly.
   opts_match <- if (requested == "cluster") {
@@ -1999,11 +2226,20 @@ pffr_vcov <- function(
   } else {
     TRUE
   }
+  adjustment_match <- requested != "cl2" ||
+    !cl2_adjustment_explicit ||
+    identical(
+      cl2_adjustment,
+      object$pffr$sandwich_info$cl2_adjustment %||%
+        object$pffr$cl2_adjustment %||%
+        "shortcut"
+    )
   if (
     !ablation &&
       is.null(cluster) &&
       identical(requested, canon$fit_type) &&
       opts_match &&
+      adjustment_match &&
       !is.null(canon$Vsandwich)
   ) {
     return(
@@ -2024,6 +2260,17 @@ pffr_vcov <- function(
   key <- if (!ablation && is.null(cluster)) {
     if (requested == "cluster" && (freq || dof_correction != "none")) {
       paste(requested, freq, dof_correction, edf_type, sep = "|")
+    } else if (requested == "cl2") {
+      if (cl2_adjustment == "auto") {
+        if (freq) paste(requested, "freq", sep = "|") else requested
+      } else {
+        paste(
+          requested,
+          cl2_adjustment,
+          if (freq) "freq" else "bayes",
+          sep = "|"
+        )
+      }
     } else if (freq) {
       paste(requested, "freq", sep = "|")
     } else {
@@ -2054,7 +2301,8 @@ pffr_vcov <- function(
     dof_correction = dof_correction,
     edf_type = edf_type,
     b2 = b2,
-    center_scores = center_scores
+    center_scores = center_scores,
+    cl2_adjustment = cl2_adjustment
   )
   if (!is.null(cache) && !is.null(key)) {
     cache[[key]] <- V
@@ -2956,6 +3204,8 @@ pffr_upgrade_fit <- function(object) {
 #'   warning) for `"cl2"`, which already corrects per-cluster leverage.
 #' @param edf_type Which EDF the `"edf"` correction uses (`"trace"`/`"edf2"`/
 #'   `"basis"`).
+#' @param cl2_adjustment CL2 leverage adjustment (`"auto"` (default),
+#'   `"exact"`, or `"shortcut"`), used only for `type = "cl2"`.
 #' @returns Model with the robust covariance stored in `$pffr$Vsandwich`.
 #' @keywords internal
 apply_sandwich_correction <- function(
@@ -2963,7 +3213,8 @@ apply_sandwich_correction <- function(
   algorithm,
   type = "cluster",
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 ) {
   gam_obj <- if (as.character(algorithm) %in% c("gamm4", "gamm")) m$gam else m
 
@@ -2992,7 +3243,8 @@ apply_sandwich_correction <- function(
     cluster_id,
     freq = FALSE,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
   Vsw_freq <- pffr_compute_sandwich(
     bread,
@@ -3000,11 +3252,14 @@ apply_sandwich_correction <- function(
     cluster_id,
     freq = TRUE,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
 
   n_capped <- attr(Vsw, "n_capped_clusters") %||% 0L
   max_lev <- attr(Vsw, "max_leverage") %||% NA_real_
+  n_adjusted <- attr(Vsw, "n_adjusted") %||% 0L
+  resolved_adjustment <- attr(Vsw, "cl2_adjustment") %||% NA_character_
 
   gam_obj$pffr$Vsandwich <- Vsw
   gam_obj$pffr$Vsandwich_freq <- Vsw_freq
@@ -3014,6 +3269,10 @@ apply_sandwich_correction <- function(
     G = if (!is.null(cluster_id)) length(unique(cluster_id)) else NA_integer_,
     n_capped = n_capped,
     max_leverage = max_lev,
+    cl2_adjustment = resolved_adjustment,
+    n_adjusted = n_adjusted,
+    min_block_eig = attr(Vsw, "min_block_eig") %||% NA_real_,
+    max_block_kappa = attr(Vsw, "max_block_kappa") %||% NA_real_,
     dof_correction = if (type == "cluster") dof_correction else "none",
     edf_type = edf_type,
     version = as.character(utils::packageVersion("refund")),
@@ -3021,6 +3280,8 @@ apply_sandwich_correction <- function(
   )
   # Keep the legacy CL2 leverage-cap diagnostic slot populated.
   gam_obj$pffr$cl2_n_capped <- if (type == "cl2") n_capped else NULL
+  gam_obj$pffr$cl2_adjustment <- if (type == "cl2") resolved_adjustment else
+    NULL
   # Fresh cache for on-demand recomputation of other sandwich types
   # (fit$pffr$Vsandwich_cache[[type]]).
   gam_obj$pffr$Vsandwich_cache <- new.env(parent = emptyenv())
diff --git a/R/pffr-methods.R b/R/pffr-methods.R
index 05736df3..3324aecc 100755
--- a/R/pffr-methods.R
+++ b/R/pffr-methods.R
@@ -1291,6 +1291,11 @@ compute_pointwise_ci <- function(
 #' @param edf_type Which EDF the \code{"edf"} correction uses
 #'   (\code{"trace"}/\code{"edf2"}/\code{"basis"}; see \code{\link{pffr}}).
 #'   Defaults to \code{NULL} (inherit from the fit).
+#' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}.
+#'   \code{NULL} (default) inherits the fit-time choice; \code{"auto"}
+#'   applies the documented feasibility rule, while \code{"exact"} and
+#'   \code{"shortcut"} force either CL2 variant. Supplying a value forces a
+#'   covariance recomputation.
 #' @param seWithMean logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in  standard errors returned for smooth components?
 #' @param n1 see below
 #' @param n2 see below
@@ -1368,6 +1373,7 @@ coef.pffr <- function(
   cluster = NULL,
   dof_correction = NULL,
   edf_type = NULL,
+  cl2_adjustment = NULL,
   seWithMean = TRUE,
   n1 = 100,
   n2 = 40,
@@ -1398,6 +1404,12 @@ coef.pffr <- function(
   if (is.null(edf_type)) edf_type <- model_edf_type
   dof_correction <- match.arg(dof_correction, c("none", "edf"))
   edf_type <- match.arg(edf_type, c("trace", "edf2", "basis"))
+  if (!is.null(cl2_adjustment)) {
+    cl2_adjustment <- match.arg(
+      cl2_adjustment,
+      c("auto", "exact", "shortcut")
+    )
+  }
   # Only warn when the user *explicitly* asked for a dof correction on a
   # non-cluster sandwich; an inherited "edf" (from the fit) stays silent.
   if (dof_explicitly_set && dof_correction != "none" && sandwich != "cluster") {
@@ -1586,7 +1598,8 @@ coef.pffr <- function(
       dof_correction = dof_correction,
       edf_type = edf_type,
       b2 = b2,
-      center_scores = center_scores
+      center_scores = center_scores,
+      cl2_adjustment = cl2_adjustment
     )
 
     # Pointwise critical-value reference (S3). `crit` selects the pointwise
diff --git a/R/pffr.R b/R/pffr.R
index 0f6668da..8f71ce73 100644
--- a/R/pffr.R
+++ b/R/pffr.R
@@ -180,7 +180,9 @@
 #'   handles both heteroskedasticity and within-curve autocorrelation — the
 #'   recommended choice for functional data when a fixed estimator is wanted.
 #'   \code{"cl2"}: leverage-adjusted cluster-robust sandwich (Bell-McCaffrey
-#'   style CL2), mainly relevant in smaller samples.
+#'   style CL2), mainly relevant in smaller samples. Within CL2, the exact
+#'   block calculation is used by default where feasible; see
+#'   \code{cl2_adjustment} below.
 #'   \code{"hc"}: observation-level HC sandwich via
 #'   \code{\link[mgcv]{vcov.gam}(sandwich = TRUE)}, which corrects for
 #'   heteroskedasticity but ignores within-curve correlation.
@@ -220,6 +222,14 @@
 #'   the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or
 #'   \code{"basis"} (the basis dimension). Only relevant when
 #'   \code{dof_correction = "edf"}.
+#' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}:
+#'   \code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
+#'   dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
+#'   10^9} (a measured sub-second marginal cost); it otherwise falls back to
+#'   the historical shortcut. The exact block
+#'   is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
+#'   \eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
+#'   variant. Ignored unless the resolved sandwich is \code{"cl2"}.
 #' @param ... additional arguments that are valid for \code{\link[mgcv]{gam}},
 #'   \code{\link[mgcv]{bam}}, \code{'\link[gamm4]{gamm4}'} or
 #'   \code{'\link[mgcv]{jagam}'}. \code{subset} is not implemented.
@@ -320,6 +330,7 @@ pffr <- function(
   sandwich = c("auto", "cluster", "cl2", "hc", "none"),
   dof_correction = c("none", "edf"),
   edf_type = c("trace", "edf2", "basis"),
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   ...
 ) {
   call <- match.call()
@@ -330,6 +341,7 @@ pffr <- function(
   sandwich <- match.arg(sandwich)
   dof_correction <- match.arg(dof_correction)
   edf_type <- match.arg(edf_type)
+  cl2_adjustment <- match.arg(cl2_adjustment)
   # "auto" is resolved after fitting (needs G / max D_g); defer the CR1 dof
   # compatibility check to then, so a legitimate auto -> "cluster" keeps its dof.
   if (dof_correction != "none" && !(sandwich %in% c("cluster", "auto"))) {
@@ -485,7 +497,8 @@ pffr <- function(
     ydata = prep$ydata,
     sandwich = sandwich,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
 
   m <- pffr_attach_metadata(m, prep$algorithm, ret)
@@ -503,6 +516,7 @@ pffr <- function(
     prep$algorithm,
     type = sandwich,
     dof_correction = dof_correction,
-    edf_type = edf_type
+    edf_type = edf_type,
+    cl2_adjustment = cl2_adjustment
   )
 }
diff --git a/man/apply_sandwich_correction.Rd b/man/apply_sandwich_correction.Rd
index 2a497f12..904a494d 100644
--- a/man/apply_sandwich_correction.Rd
+++ b/man/apply_sandwich_correction.Rd
@@ -9,7 +9,8 @@ apply_sandwich_correction(
   algorithm,
   type = "cluster",
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 )
 }
 \arguments{
@@ -26,6 +27,9 @@ warning) for `"cl2"`, which already corrects per-cluster leverage.}
 
 \item{edf_type}{Which EDF the `"edf"` correction uses (`"trace"`/`"edf2"`/
 `"basis"`).}
+
+\item{cl2_adjustment}{CL2 leverage adjustment (`"auto"` (default),
+`"exact"`, or `"shortcut"`), used only for `type = "cl2"`.}
 }
 \value{
 Model with the robust covariance stored in `$pffr$Vsandwich`.
diff --git a/man/coef.pffr.Rd b/man/coef.pffr.Rd
index 11eceb3b..a6e8d4de 100644
--- a/man/coef.pffr.Rd
+++ b/man/coef.pffr.Rd
@@ -13,6 +13,7 @@
   cluster = NULL,
   dof_correction = NULL,
   edf_type = NULL,
+  cl2_adjustment = NULL,
   seWithMean = TRUE,
   n1 = 100,
   n2 = 40,
@@ -68,6 +69,12 @@ Ignored (with a warning) for \code{sandwich} other than \code{"cluster"}.}
 (\code{"trace"}/\code{"edf2"}/\code{"basis"}; see \code{\link{pffr}}).
 Defaults to \code{NULL} (inherit from the fit).}
 
+\item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}.
+\code{NULL} (default) inherits the fit-time choice; \code{"auto"}
+applies the documented feasibility rule, while \code{"exact"} and
+\code{"shortcut"} force either CL2 variant. Supplying a value forces a
+covariance recomputation.}
+
 \item{seWithMean}{logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in  standard errors returned for smooth components?}
 
 \item{n1}{see below}
diff --git a/man/gam_sandwich_cluster_cl2.Rd b/man/gam_sandwich_cluster_cl2.Rd
index f5d1ac61..ee0668f0 100644
--- a/man/gam_sandwich_cluster_cl2.Rd
+++ b/man/gam_sandwich_cluster_cl2.Rd
@@ -10,6 +10,7 @@ gam_sandwich_cluster_cl2(
   freq = FALSE,
   tol = 1e-08,
   leverage_cap = 0.999,
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   b2 = TRUE,
   center_scores = FALSE
 )
@@ -27,6 +28,11 @@ If `FALSE` (default), use Bayesian sandwich (`B2 = Vp - Ve`).}
 
 \item{leverage_cap}{Cap for cluster leverage eigenvalues (< 1).}
 
+\item{cl2_adjustment}{CL2 leverage adjustment: `"auto"` (default) selects
+the exact Bell--McCaffrey block where it is relevant and affordable,
+`"exact"` forces that block, and `"shortcut"` uses the historical
+per-cluster shortcut `(I - H_gg)^(-1/2)`.}
+
 \item{b2}{Internal ablation switch (default `TRUE` = current behavior). When
 `FALSE`, drop the additive Bayesian smoothing-bias term \eqn{B_2 = V_p -
 V_e} (X5).}
@@ -37,8 +43,13 @@ contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat
 (X6).}
 }
 \value{
-A p x p covariance matrix with attributes `n_capped_clusters` and
-  `max_leverage` for the CL2 leverage diagnostic.
+A p x p covariance matrix with leverage diagnostics. Exact CL2
+  returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`),
+  `min_block_eig`, `max_block_kappa` and `min_block_eig_rel`; the shortcut
+  returns the legacy `n_capped_clusters` and `max_leverage` attributes. Both
+  return `max_obs_leverage` and `hat_invariant_violation` (`NULL` when the
+  penalized hat respects its bounds, otherwise a description of the
+  violation; see [pffr_hat_invariant_violation()]).
 }
 \description{
 Computes a cluster-robust covariance matrix with a Bell-McCaffrey style
diff --git a/man/pffr.Rd b/man/pffr.Rd
index 4fe62cae..1521b53f 100644
--- a/man/pffr.Rd
+++ b/man/pffr.Rd
@@ -17,6 +17,7 @@ pffr(
   sandwich = c("auto", "cluster", "cl2", "hc", "none"),
   dof_correction = c("none", "edf"),
   edf_type = c("trace", "edf2", "basis"),
+  cl2_adjustment = c("auto", "exact", "shortcut"),
   ...
 )
 }
@@ -81,7 +82,9 @@ penalty.}
   handles both heteroskedasticity and within-curve autocorrelation — the
   recommended choice for functional data when a fixed estimator is wanted.
   \code{"cl2"}: leverage-adjusted cluster-robust sandwich (Bell-McCaffrey
-  style CL2), mainly relevant in smaller samples.
+  style CL2), mainly relevant in smaller samples. Within CL2, the exact
+  block calculation is used by default where feasible; see
+  \code{cl2_adjustment} below.
   \code{"hc"}: observation-level HC sandwich via
   \code{\link[mgcv]{vcov.gam}(sandwich = TRUE)}, which corrects for
   heteroskedasticity but ignores within-curve correlation.
@@ -124,6 +127,15 @@ the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or
 \code{"basis"} (the basis dimension). Only relevant when
 \code{dof_correction = "edf"}.}
 
+\item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}:
+\code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the
+dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times
+10^9} (a measured sub-second marginal cost); it otherwise falls back to
+the historical shortcut. The exact block
+is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at
+\eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either
+variant. Ignored unless the resolved sandwich is \code{"cl2"}.}
+
 \item{...}{additional arguments that are valid for \code{\link[mgcv]{gam}},
 \code{\link[mgcv]{bam}}, \code{'\link[gamm4]{gamm4}'} or
 \code{'\link[mgcv]{jagam}'}. \code{subset} is not implemented.}
diff --git a/man/pffr_build_metadata.Rd b/man/pffr_build_metadata.Rd
index 925ec667..9d183451 100644
--- a/man/pffr_build_metadata.Rd
+++ b/man/pffr_build_metadata.Rd
@@ -24,7 +24,8 @@ pffr_build_metadata(
   ydata,
   sandwich,
   dof_correction = "none",
-  edf_type = "trace"
+  edf_type = "trace",
+  cl2_adjustment = "auto"
 )
 }
 \description{
diff --git a/man/pffr_compute_sandwich.Rd b/man/pffr_compute_sandwich.Rd
index 5d60df7e..6069f230 100644
--- a/man/pffr_compute_sandwich.Rd
+++ b/man/pffr_compute_sandwich.Rd
@@ -12,7 +12,8 @@ pffr_compute_sandwich(
   dof_correction = "none",
   edf_type = "trace",
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = "auto"
 )
 }
 \arguments{
@@ -34,6 +35,9 @@ only).}
 \item{center_scores}{Internal ablation switch (default `FALSE`); center the
 per-cluster score sums before the meat when `TRUE` (X6). Applies to
 `"cluster"`/`"cl2"` only.}
+
+\item{cl2_adjustment}{CL2 leverage adjustment (`"auto"`, `"exact"`, or
+`"shortcut"`), used only for `type = "cl2"`.}
 }
 \value{
 A covariance matrix (with CL2 leverage attributes for `type =
diff --git a/man/pffr_hat_invariant_violation.Rd b/man/pffr_hat_invariant_violation.Rd
new file mode 100644
index 00000000..3c0ed713
--- /dev/null
+++ b/man/pffr_hat_invariant_violation.Rd
@@ -0,0 +1,54 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pffr-core.R
+\name{pffr_hat_invariant_violation}
+\alias{pffr_hat_invariant_violation}
+\title{Numerical-sanity check on the per-cluster hat blocks}
+\usage{
+pffr_hat_invariant_violation(
+  max_obs_leverage = NA_real_,
+  max_leverage = NA_real_,
+  min_block_eig_rel = 0,
+  tol = 1e-06
+)
+}
+\arguments{
+\item{max_obs_leverage}{Largest per-observation leverage \eqn{h_{ii}} seen,
+or `NA`.}
+
+\item{max_leverage}{Largest eigenvalue of any \eqn{H_{gg}} (shortcut path
+only; `NA` otherwise).}
+
+\item{min_block_eig_rel}{Smallest eigenvalue of any \eqn{B_g} relative to
+that block's largest eigenvalue (exact path only; `0` otherwise).}
+
+\item{tol}{Relative slack allowed before an invariant counts as violated.}
+}
+\value{
+`NULL` when every invariant holds, otherwise a one-sentence
+  character description of the violation.
+}
+\description{
+The penalized hat \eqn{H = X_w V_p X_w'} obeys \eqn{0 \le h_{ii} \le 1} and,
+since any principal block of a symmetric matrix has its eigenvalues inside
+the parent's range, \eqn{0 \le \mathrm{eigen}(H_{gg}) \le 1}. The exact
+Bell--McCaffrey block \eqn{B_g = ((I - H)^2)_{gg}} is positive semi-definite
+for the same reason. When a fit is numerically degenerate --- an
+ill-conditioned or barely converged fit, e.g. a Poisson fit whose fitted
+means span many orders of magnitude --- the bread \eqn{V_p} stops being the
+inverse of the same weighted cross-product, and these invariants break by
+orders of magnitude rather than by round-off. This was the mechanism behind
+the numerically exploded interval widths observed on hard Poisson
+replicates in the locked benchmark (study LB, claim P-LB5).
+}
+\details{
+Detection is deliberately separate from the leverage cap. Capping
+\eqn{H_{gg}} at `leverage_cap` silently converts an impossible leverage into
+the *largest legitimate* one, so the shortcut adjustment responds to garbage
+with its maximal variance inflation; the exact block, being a squared
+operator, instead deflates. Both are bounded, neither is meaningful, and
+only a diagnostic can tell the user which situation they are in. We
+therefore warn and keep going rather than substituting a differently wrong
+number: the CR1 covariance is built from the same bread and is no more
+trustworthy here.
+}
+\keyword{internal}
diff --git a/man/pffr_vcov.Rd b/man/pffr_vcov.Rd
index 051bc665..697cbe58 100644
--- a/man/pffr_vcov.Rd
+++ b/man/pffr_vcov.Rd
@@ -12,7 +12,8 @@ pffr_vcov(
   dof_correction = NULL,
   edf_type = NULL,
   b2 = TRUE,
-  center_scores = FALSE
+  center_scores = FALSE,
+  cl2_adjustment = NULL
 )
 }
 \arguments{
@@ -39,6 +40,9 @@ recomputation and are never served from or written to the cache.}
 \item{center_scores}{Internal ablation switch (default `FALSE` = current
 behavior). When `TRUE`, center the per-cluster score sums before forming
 the cluster/CL2 meat (X6). Same cache semantics as `b2`.}
+
+\item{cl2_adjustment}{CL2 leverage adjustment. `NULL` (default) inherits
+the fitted choice; otherwise one of `"auto"`, `"exact"`, or `"shortcut"`.}
 }
 \value{
 A covariance matrix.
diff --git a/man/resolve_cl2_adjustment.Rd b/man/resolve_cl2_adjustment.Rd
new file mode 100644
index 00000000..009b4c2f
--- /dev/null
+++ b/man/resolve_cl2_adjustment.Rd
@@ -0,0 +1,40 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pffr-core.R
+\name{resolve_cl2_adjustment}
+\alias{resolve_cl2_adjustment}
+\title{Resolve the CL2 leverage adjustment}
+\usage{
+resolve_cl2_adjustment(
+  cl2_adjustment = c("auto", "exact", "shortcut"),
+  G,
+  maxDg,
+  p
+)
+}
+\arguments{
+\item{cl2_adjustment}{One of `"auto"`, `"exact"`, or `"shortcut"`.}
+
+\item{G}{Number of clusters.}
+
+\item{maxDg}{Largest cluster size.}
+
+\item{p}{Coefficient-space dimension.}
+}
+\value{
+`"exact"` or `"shortcut"`.
+}
+\description{
+Exact CL2 is the default when the finite-sample correction can matter
+(`G <= 100`; Study EX measured the exact-over-shortcut coverage gain
+decaying from ~+1pp at G = 20 to under +0.15pp at G = 100) and the dense
+block multiplication is affordable. The cost proxy `G * maxDg * p^2`
+covers the dominant `(Xw_g Vp) (Xw'Xw)` multiplication. The 5e9 bound is
+calibrated by measurement, not guessed: both adjustments pay the same
+per-cluster eigendecomposition, so the exact block's MARGINAL cost is only
+these multiplications -- timed at ~0.03 s extra at proxy 5e8 and ~0.7 s
+total at 5e9 (single-thread BLAS), i.e. the bound caps the correctness
+upgrade at well under a second and in particular never denies it in the
+saturated small-G regime where its benefit is largest. Large-G CL2 uses
+the historical shortcut because the correction tends to zero there.
+}
+\keyword{internal}
diff --git a/tests/testthat/helper-exactcl2.R b/tests/testthat/helper-exactcl2.R
new file mode 100644
index 00000000..711e7bd2
--- /dev/null
+++ b/tests/testthat/helper-exactcl2.R
@@ -0,0 +1,158 @@
+# Fixtures for the exact-CL2 hardening tests (DESIGN amendments A18 and the
+# 2026-08-04 study-LB P-LB5 follow-up). The full standalone gates, which also
+# write the archived result tables, live in the pffr-ci repo
+# (diagnostics/hardening-exactcl2.R); these helpers keep the same coverage
+# available to CI without a non-standard top-level file in the package.
+
+make_exactcl2_fixture <- function(family_name, k, design) {
+  set.seed(
+    match(family_name, c("gaussian", "poisson", "binomial", "Gamma")) *
+      100 +
+      k +
+      match(design, c("unbalanced", "influential"))
+  )
+  n_curve <- 12L
+  n_grid <- 10L
+  grid <- seq(0, 1, length.out = n_grid)
+  xlin <- stats::rnorm(n_curve)
+  eta <- outer(xlin, rep(0.35, n_grid)) +
+    matrix(rep(0.2 * sin(2 * pi * grid), each = n_curve), n_curve, n_grid) -
+    0.4
+  Y <- switch(
+    family_name,
+    gaussian = eta + matrix(stats::rnorm(n_curve * n_grid, sd = 0.35), n_curve),
+    poisson = matrix(
+      stats::rpois(n_curve * n_grid, lambda = exp(pmin(eta, 2))),
+      n_curve
+    ),
+    binomial = matrix(
+      stats::rbinom(n_curve * n_grid, 1, stats::plogis(eta)),
+      n_curve
+    ),
+    Gamma = matrix(
+      stats::rgamma(n_curve * n_grid, shape = 5, scale = exp(pmin(eta, 2)) / 5),
+      n_curve
+    )
+  )
+  # Keep the response in its ordinary range, then make one covariate curve
+  # extreme. This isolates high leverage from a separable/outcome-saturation
+  # failure and exercises the shipped eigenvalue floor.
+  if (design == "influential") xlin[1] <- 3000
+  # Five deliberately unbalanced independent clusters, with cluster 1 made
+  # influential above. Their curve counts 1, 1, 2, 3, 5 yield unequal blocks.
+  cluster <- rep(seq_len(5), times = c(1, 1, 2, 3, 5))
+  list(
+    data = data.frame(Y = I(Y), xlin = xlin),
+    yind = grid,
+    cluster = cluster,
+    family = switch(
+      family_name,
+      gaussian = stats::gaussian(),
+      poisson = stats::poisson(),
+      binomial = stats::binomial(),
+      Gamma = stats::Gamma(link = "log")
+    ),
+    k = k
+  )
+}
+
+run_exactcl2_hardening <- function() {
+  cases <- expand.grid(
+    family = c("gaussian", "poisson", "binomial", "Gamma"),
+    k = c(4L, 7L),
+    design = c("unbalanced", "influential"),
+    stringsAsFactors = FALSE
+  )
+  rows <- vector("list", nrow(cases))
+
+  for (i in seq_len(nrow(cases))) {
+    case <- cases[i, ]
+    fixture <- make_exactcl2_fixture(case$family, case$k, case$design)
+    fit <- suppressWarnings(suppressMessages(refund::pffr(
+      Y ~ xlin,
+      data = fixture$data,
+      yind = fixture$yind,
+      family = fixture$family,
+      bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+      sandwich = "none"
+    )))
+
+    # `1 - sqrt(tol)` makes the exact floor equal tol: Study EX's
+    # cl2_exact_nocap comparator without duplicating the implementation.
+    V_exact <- suppressWarnings(refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cluster = fixture$cluster,
+      cl2_adjustment = "exact"
+    ))
+    V_nocap <- suppressWarnings(refund:::gam_sandwich_cluster_cl2(
+      refund:::pffr_model_based_gam(fit),
+      refund:::build_cluster_id(fit$pffr, cluster = fixture$cluster),
+      cl2_adjustment = "exact",
+      leverage_cap = 1 - sqrt(1e-8)
+    ))
+    V_shortcut <- suppressWarnings(refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cluster = fixture$cluster,
+      cl2_adjustment = "shortcut"
+    ))
+
+    se_exact <- sqrt(pmax(diag(V_exact), 0))
+    se_shortcut <- sqrt(pmax(diag(V_shortcut), 0))
+    positive <- is.finite(se_exact) &
+      is.finite(se_shortcut) &
+      se_exact > 0 &
+      se_shortcut > 0
+    ratio <- se_exact[positive] / se_shortcut[positive]
+    rows[[i]] <- data.frame(
+      family = case$family,
+      k = case$k,
+      design = case$design,
+      G = length(unique(fixture$cluster)),
+      max_cluster_size = max(table(fixture$cluster)) * length(fixture$yind),
+      n_adjusted = attr(V_exact, "n_adjusted"),
+      min_block_eig = attr(V_exact, "min_block_eig"),
+      max_block_kappa = attr(V_exact, "max_block_kappa"),
+      max_abs_exact_nocap = max(abs(V_exact - V_nocap)),
+      median_se_ratio_to_shortcut = stats::median(ratio),
+      se_finite_positive = all(positive),
+      se_ratio_sane = length(ratio) > 0 &&
+        min(ratio) > 1e-3 &&
+        max(ratio) < 1e3,
+      stringsAsFactors = FALSE
+    )
+  }
+  do.call(rbind, rows)
+}
+
+#' Study-LB P-LB5 fixture: a Poisson design whose latent roughness amplitude
+#' `amp` pushes individual fitted means far above the nominal marginal mean of
+#' 4, the mechanism behind the exploded interval widths in the locked
+#' benchmark. `amp = 1` is the benign control; `amp = 10` is degenerate.
+make_lb5_fixture <- function(amp, n_grid, k, seed) {
+  set.seed(seed)
+  n_curve <- 8L
+  grid <- seq(0, 1, length.out = n_grid)
+  xlin <- stats::rnorm(n_curve)
+  nfpc <- 6L
+  rough <- matrix(stats::rnorm(n_curve * nfpc), n_curve, nfpc) %*%
+    t(sapply(seq_len(nfpc), function(j) sin(j * pi * grid)))
+  eta <- outer(xlin, rep(0.35, n_grid)) + amp * rough + log(4)
+  Y <- matrix(
+    stats::rpois(n_curve * n_grid, lambda = exp(pmin(eta, 600))),
+    n_curve
+  )
+  list(data = data.frame(Y = I(Y), xlin = xlin), yind = grid, k = k)
+}
+
+fit_lb5_fixture <- function(fixture) {
+  suppressWarnings(suppressMessages(refund::pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    family = stats::poisson(),
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "none"
+  )))
+}
diff --git a/tests/testthat/test-pffr-exactcl2.R b/tests/testthat/test-pffr-exactcl2.R
new file mode 100644
index 00000000..d9b88b10
--- /dev/null
+++ b/tests/testthat/test-pffr-exactcl2.R
@@ -0,0 +1,143 @@
+# A18 exact-CL2 hardening gate. The standalone script writes the same table
+# for release records; this fast test keeps its family/design coverage in CI.
+# Fixtures live in helper-exactcl2.R (they used to be sourced from a top-level
+# script that was removed from the package tree, which left this file unable to
+# run at all).
+
+test_that("exact CL2 is robust on high-leverage family and basis fixtures", {
+  results <- run_exactcl2_hardening()
+
+  expect_equal(nrow(results), 16L)
+  expect_true(all(is.finite(results$n_adjusted)))
+  expect_gt(sum(results$n_adjusted), 0)
+  expect_true(all(is.finite(results$min_block_eig)))
+  expect_true(all(is.finite(results$max_block_kappa)))
+  expect_true(all(is.finite(results$max_abs_exact_nocap)))
+  expect_true(all(results$se_finite_positive))
+  expect_true(all(results$se_ratio_sane))
+})
+
+test_that("CL2 adjustment selection can be forced and is recorded", {
+  fixture <- make_exactcl2_fixture("gaussian", 4L, "influential")
+  fit_auto <- suppressWarnings(suppressMessages(pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "cl2"
+  )))
+  fit_exact <- suppressWarnings(suppressMessages(pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "cl2",
+    cl2_adjustment = "exact"
+  )))
+  fit_shortcut <- suppressWarnings(suppressMessages(pffr(
+    Y ~ xlin,
+    data = fixture$data,
+    yind = fixture$yind,
+    bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)),
+    sandwich = "cl2",
+    cl2_adjustment = "shortcut"
+  )))
+
+  expect_identical(fit_exact$pffr$sandwich_info$cl2_adjustment, "exact")
+  expect_identical(fit_shortcut$pffr$sandwich_info$cl2_adjustment, "shortcut")
+  expect_identical(fit_auto$pffr$sandwich_info$cl2_adjustment, "exact")
+  expect_equal(fit_auto$pffr$Vsandwich, fit_exact$pffr$Vsandwich)
+  expect_true(is.numeric(fit_exact$pffr$sandwich_info$n_adjusted))
+  expect_true(is.numeric(fit_shortcut$pffr$sandwich_info$n_adjusted))
+})
+
+test_that("automatic exact CL2 uses the documented relevance and cost rule", {
+  resolve <- refund:::resolve_cl2_adjustment
+  expect_identical(resolve("auto", G = 100, maxDg = 50, p = 20), "exact")
+  expect_identical(resolve("auto", G = 101, maxDg = 50, p = 20), "shortcut")
+  # 9e8 ops: below the 5e9 cost cap (raised from 5e8, see notes 2026-07-21)
+  expect_identical(resolve("auto", G = 20, maxDg = 500, p = 300), "exact")
+  # 6.4e10 ops: above the cap even at eligible G
+  expect_identical(resolve("auto", G = 100, maxDg = 1000, p = 800), "shortcut")
+  expect_identical(resolve("exact", G = 1000, maxDg = 5000, p = 500), "exact")
+  expect_identical(resolve("shortcut", G = 2, maxDg = 2, p = 2), "shortcut")
+})
+
+# ---------------------------------------------------------------------------
+# Study-LB P-LB5 follow-up (2026-08-04): numerically exploded interval widths
+# on degenerate Poisson fits must be DETECTED, not returned silently.
+# ---------------------------------------------------------------------------
+
+test_that("pffr_hat_invariant_violation() flags only impossible hat values", {
+  viol <- refund:::pffr_hat_invariant_violation
+
+  # Every invariant holds -> NULL.
+  expect_null(viol(max_obs_leverage = 0.99, max_leverage = 0.999))
+  expect_null(viol(max_obs_leverage = NA_real_, max_leverage = NA_real_))
+  # Round-off-sized excess is tolerated.
+  expect_null(viol(max_obs_leverage = 1 + 1e-9))
+  # Real violations are described, with the offending number in the message.
+  expect_match(
+    viol(max_obs_leverage = 22.6),
+    "per-observation leverage is 22.6"
+  )
+  expect_match(viol(max_leverage = 8.18e84), "hat eigenvalue")
+  expect_match(viol(min_block_eig_rel = -0.3), "positive semi-definite")
+  # A non-finite monitor must not trip the check.
+  expect_null(viol(max_obs_leverage = Inf, max_leverage = NaN))
+})
+
+test_that("a benign Poisson CL2 fit reports leverages inside the bounds", {
+  skip_on_cran()
+  fit <- fit_lb5_fixture(make_lb5_fixture(amp = 1, n_grid = 20L, k = 8L, 21L))
+  V_exact <- expect_no_warning(
+    refund:::pffr_vcov(fit, sandwich = "cl2", cl2_adjustment = "exact")
+  )
+  V_short <- expect_no_warning(
+    refund:::pffr_vcov(fit, sandwich = "cl2", cl2_adjustment = "shortcut")
+  )
+  expect_lte(attr(V_exact, "max_obs_leverage"), 1)
+  expect_lte(attr(V_short, "max_leverage"), 1)
+  expect_null(attr(V_exact, "hat_invariant_violation"))
+  expect_null(attr(V_short, "hat_invariant_violation"))
+  # exact and shortcut agree closely in the well-behaved regime
+  ratio <- sqrt(max(abs(diag(V_exact)))) / sqrt(max(abs(diag(V_short))))
+  expect_gt(ratio, 0.5)
+  expect_lt(ratio, 2)
+})
+
+test_that("a degenerate Poisson fit warns instead of returning a silent 1e30", {
+  skip_on_cran()
+  # amp = 10 drives individual fitted means many orders of magnitude above the
+  # nominal marginal mean, the study-LB mechanism; the penalized hat then
+  # breaks its own bounds by dozens of orders of magnitude.
+  fit <- fit_lb5_fixture(make_lb5_fixture(amp = 10, n_grid = 30L, k = 12L, 21L))
+
+  expect_warning(
+    V_exact <- refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cl2_adjustment = "exact"
+    ),
+    "NOT trustworthy"
+  )
+  expect_warning(
+    V_short <- refund:::pffr_vcov(
+      fit,
+      sandwich = "cl2",
+      cl2_adjustment = "shortcut"
+    ),
+    "NOT trustworthy"
+  )
+  expect_gt(attr(V_exact, "max_obs_leverage"), 1)
+  expect_type(attr(V_exact, "hat_invariant_violation"), "character")
+
+  se_exact <- sqrt(max(abs(diag(V_exact))))
+  se_short <- sqrt(max(abs(diag(V_short))))
+  expect_true(is.finite(se_exact) && is.finite(se_short))
+  # The exact block is a principal block of a squared symmetric matrix and so
+  # stays positive semi-definite; the shortcut's I - H_gg turns indefinite and
+  # is floored into the largest legitimate variance inflation. The exact path
+  # is therefore far less explosive here, though neither is meaningful.
+  expect_lt(se_exact, se_short)
+})
diff --git a/tests/testthat/test-pffr.R b/tests/testthat/test-pffr.R
index c9dca72c..37ad8d3a 100644
--- a/tests/testthat/test-pffr.R
+++ b/tests/testthat/test-pffr.R
@@ -1855,17 +1855,18 @@ test_that("CL2 reports leverage diagnostics and warns when the cap is hit", {
     V_benign <- gam_sandwich_cluster_cl2(b, cluster_id, freq = FALSE),
     NA
   )
-  expect_identical(attr(V_benign, "n_capped_clusters"), 0L)
-  expect_true(is.finite(attr(V_benign, "max_leverage")))
-  expect_gt(attr(V_benign, "max_leverage"), 0)
-  expect_lt(attr(V_benign, "max_leverage"), 0.999)
+  expect_identical(attr(V_benign, "cl2_adjustment"), "exact")
+  expect_identical(attr(V_benign, "n_adjusted"), 0L)
+  expect_true(is.finite(attr(V_benign, "min_block_eig")))
+  expect_gt(attr(V_benign, "min_block_eig"), 0)
 
   expect_warning(
     V_capped <- gam_sandwich_cluster_cl2(
       b,
       cluster_id,
       freq = FALSE,
-      leverage_cap = 0.001
+      leverage_cap = 0.001,
+      cl2_adjustment = "shortcut"
     ),
     "CL2 leverage adjustment hit the leverage cap"
   )
@@ -1876,7 +1877,8 @@ test_that("CL2 reports leverage diagnostics and warns when the cap is hit", {
     b,
     cluster_id,
     freq = FALSE,
-    leverage_cap = 0.001
+    leverage_cap = 0.001,
+    cl2_adjustment = "shortcut"
   ))
   expect_equal(
     as.matrix(V_capped),