diff --git a/DESCRIPTION b/DESCRIPTION
index adc319bd8..feabaf90b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -55,10 +55,10 @@ Imports:
   utils,
 Suggests:
   cluster,
+  Coreset,
   future,
   highs,
   knitr,
-  Coreset,
   phangorn (>= 2.2.1),
   PlotTools,
   promises,
diff --git a/NEWS.md b/NEWS.md
index fb7534252..f69a23ba1 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,36 @@
 # To integrate into 2.0.0 notes
 
+- `QuartetConcordance()` gains `unit`, selecting the currency in which quartets
+  are counted.  `unit = "nrqs"` (the new default) counts only non-redundant
+  quartet statements: of the quartets resolved by a split of sizes
+  (_k_, _t_ − _k_), just (_k_ − 1)(_t_ − _k_ − 1) are non-redundant,
+  so a character now scores full marks only where it
+  displays the split, rather than wherever it agrees with a large combinatorial
+  volume of quartets.  Pass `unit = "quartet"` to recover the raw count.
+
+  Note that "full marks where it displays the split" bounds each split
+  *individually*, and so applies to `return = "edge"`.  It does **not** apply to
+  `return = "char"`, which averages over every split in the tree: a character
+  agrees exactly with at most one split and is merely compatible with -- silent
+  about -- the rest, and coverage scores silence as a failure to cover.  A
+  character identical to one of the tree's own splits therefore scores well
+  below 1 (0.49 for a 22|26 split of a 48-leaf tree, 0.25 for a 2|46 split), and
+  the attainable maximum depends on the tree.  Per-character `unit = "nrqs"`
+  values rank characters within a tree but must not be read against a ceiling of
+  1 or compared across trees; use `unit = "quartet"` where a per-character index
+  bounded at 1 is wanted.  Documentation only -- no numeric behaviour changed.
+
+- The `normalize` argument of `ClusteringConcordance()`, `ConcordanceTable()`
+  and `QuartetConcordance()` is renamed `chanceCorrect`, which is what it
+  controls: the zero point is moved to the value expected under a fixed-marginal
+  null.  (The scaling to a maximum of 1 was never governed by this argument.)
+
+- `QuartetConcordance()` now applies that chance correction by default
+  (`chanceCorrect = TRUE`), matching `ClusteringConcordance()`.  Values may
+  therefore be negative, denoting agreement below chance expectation; pass
+  `chanceCorrect = FALSE` for the uncorrected measure reported by earlier
+  versions.
+
 - New `SearchControl()` parameter `enumMaxTrees`: a retention ceiling applied to
   the post-search MPT-enumeration phase alone.  `0` (the default)
   keeps `poolMaxSize` throughout, so behaviour is unchanged unless you set it.
diff --git a/R/Concordance.R b/R/Concordance.R
index 9cb70455d..a532443a0 100644
--- a/R/Concordance.R
+++ b/R/Concordance.R
@@ -60,22 +60,6 @@ NULL
 #' so that 1 corresponds to the maximum possible mutual information for each
 #' split–character pair (`hBest`).
 #'
-#' The `normalize` argument specifies how the zero point is defined.
-#'
-#' - If `normalize = FALSE`, zero corresponds to *zero* MI, without correcting
-#'   for the positive bias that arises because MI is rarely exactly zero in
-#'   finite samples.
-#'
-#' - If `normalize = TRUE`, the expected MI is computed using an analytical
-#'   approximation based on the distribution of character tokens. This is fast
-#'   and generally accurate for large trees (~200+ taxa), but does not account 
-#'   for correlation between splits.
-#'
-#' - If `normalize` is a positive integer `n`, the expected MI is estimated
-#'   empirically by fitting each character to `n` uniformly random trees and
-#'   averaging the resulting MI values. This Monte Carlo approach provides a
-#'   more accurate baseline for small trees, for which the analytical
-#'   approximation is biased. Monte Carlo standard errors are returned.
 #'
 #' @param return Character specifying the summary to return. Options are:
 #'   - `"edge"`: average concordance of each tree split across all characters;
@@ -84,18 +68,19 @@ NULL
 #'   - `"tree"`: an overall tree‑level concordance score;
 #'   - `"all"`: a full array of MI components and normalized values for every
 #'     split–character pair.
-#' 
-#'   Matching is case‑insensitive and partial.
 #'
-#' @param normalize Controls how the *expected* mutual information (the zero
-#' point of the scale) is determined.  
-#'   - `FALSE`: no chance correction; MI is scaled only by its maximum.  
-#'   - `TRUE`: subtract the analytical expected MI for random association.  
-#'   - `<integer>`: subtract an empirical expected MI estimated from that
-#'     number of random trees.
-#' 
-#'   In all cases, 1 corresponds to the maximal attainable MI for the pair
-#'   (`hBest`), and 0 corresponds to the chosen expectation.
+#' @param chanceCorrect Sets the zero point of the scale.
+#' - If `FALSE`, zero corresponds to zero MI.
+#' - If `TRUE`, zero is the value expected when each character's tokens are
+#' reassigned at random across the leaves, holding its state frequencies and the
+#' split sizes fixed: `QuartetConcordance()` computes this expectation exactly,
+#' whereas `ClusteringConcordance()` approximates it, accurately for large trees
+#' (~200+ taxa) but neglecting correlation between splits.
+#' - If a positive integer `n`, the expectation is sampled over `n` random
+#' reassignments of each character's tokens (`QuartetConcordance()`), or against
+#' `n` uniformly random trees (`ClusteringConcordance()`).
+#' _Hint: Clamp chance-corrected values to \eqn{[-1, 1]} before plotting with [QCol()] /
+#' [QACol()]_.
 #' 
 #' @returns
 #' `ClusteringConcordance(return = "all")` returns a 3D array where each
@@ -116,7 +101,7 @@ NULL
 #' corresponds to a split (an edge of the tree) and gives the normalized mutual
 #' information between that split and the character data, averaged across all
 #' characters.
-#' When `normalize = TRUE` (default), values are scaled relative to random
+#' When `chanceCorrect = TRUE` (default), values are scaled relative to random
 #' expectation; when `FALSE`, raw mutual information normalized by `hBest` is
 #' returned.
 #'
@@ -158,7 +143,7 @@ ClusteringConcordance <- function(
   tree,
   dataset,
   return = "edge",
-  normalize = TRUE
+  chanceCorrect = TRUE
 ) {
   # Check inputs
   if (is.null(dataset)) {
@@ -255,7 +240,7 @@ ClusteringConcordance <- function(
                        hh["hSplit", , , drop = FALSE] -
                        hh["hJoint", , , drop = FALSE], NULL)
   miRand <- `rownames<-`(hh["miRand", , , drop = FALSE], NULL)
-  norm <- if (isFALSE(normalize)) {
+  norm <- if (isFALSE(chanceCorrect)) {
     ifelse(hBest == 0, NA, mi / hBest)
   } else {
     ifelse(hBest == 0, NA, .Rezero(mi / hBest, miRand / hBest))
@@ -269,8 +254,8 @@ ClusteringConcordance <- function(
     charMax <- vapply(charSplits, ClusteringEntropy, double(1))[
       attr(dataset, "index")]
     charInfo <- MutualClusteringInfo(tree, charSplits)[at[["index"]]]
-    if (is.numeric(normalize)) {
-      rTrees <- replicate(normalize, RandomTree(tree), simplify = FALSE)
+    if (is.numeric(chanceCorrect)) {
+      rTrees <- replicate(chanceCorrect, RandomTree(tree), simplify = FALSE)
       # Score each random tree against `charSplits` separately: characters with
       # ambiguous tokens yield splits over different tip subsets, and the
       # vectorised `MutualClusteringInfo(<list of trees>, <list of splits>)`
@@ -283,12 +268,12 @@ ClusteringConcordance <- function(
         double(length(charSplits))
       ))[, attr(dataset, "index"), drop = FALSE]
       randMean <- colMeans(randInfo)
-      var <- rowSums((t(randInfo) - randMean) ^ 2) / (normalize - 1)
-      mcse <- sqrt(var / normalize)
+      var <- rowSums((t(randInfo) - randMean) ^ 2) / (chanceCorrect - 1)
+      mcse <- sqrt(var / chanceCorrect)
       randTreeInfo <- rowSums(randInfo)
       randTreeMean <- mean(randTreeInfo)
       treeVar <- var(randTreeInfo)
-      mcseTree <- sqrt(treeVar / normalize)
+      mcseTree <- sqrt(treeVar / chanceCorrect)
     }
   }
 
@@ -307,7 +292,7 @@ ClusteringConcordance <- function(
            best <- rowSums(hBest[1, , , drop = FALSE], dims = 2)
            ifelse(!is.na(best) & best == 0,
                   NA_real_,
-                  if (isTRUE(normalize)) {
+                  if (isTRUE(chanceCorrect)) {
                     .Rezero(
                       rowSums(mi[1, , , drop = FALSE], dims = 2) / best,
                       rowSums(miRand[1, , , drop = FALSE], dims = 2) / best
@@ -320,27 +305,28 @@ ClusteringConcordance <- function(
          
            # one <- hh["hChar", 1, , drop = TRUE] # All rows equal
            one <- charMax
-           zero <- if (isFALSE(normalize)) {
+           zero <- if (isFALSE(chanceCorrect)) {
              0
-           } else if (isTRUE(normalize)) {
+           } else if (isTRUE(chanceCorrect)) {
              apply(.ConcSlice(hh, "miRand"), 2, max)
            } else {
              randMean
            }
            ret <- (charInfo - zero) / (one - zero)
-           if (is.numeric(normalize)) {
+           if (is.numeric(chanceCorrect)) {
              mcseInfo <- ((one - charInfo) / (one - zero) ^ 2) * mcse
              mcseInfo[mcseInfo < sqrt(.Machine$double.eps)] <- 0
              structure(ret, hMax = charMax, mcse = mcseInfo)
            } else {
              # The characterwise return is deliberately NOT random-expectation
-             # normalized for logical `normalize`: `charInfo` is
+             # normalized for logical `chanceCorrect`: `charInfo` is
              # MutualClusteringInfo() against the whole tree, whereas the
              # analytic `zero` baseline above is per-single-split expected MI, so
              # subtracting it would mix incompatible quantities (and the
              # entropy-weighted variant was abandoned -- see the note below the
-             # @return docs). Only the Monte-Carlo path (numeric `normalize`)
-             # offers a same-scale empirical baseline. So return charInfo scaled
+             # @return docs). Only the Monte-Carlo path (numeric
+             # `chanceCorrect`) offers a same-scale empirical baseline, so we
+             # return charInfo scaled
              # by its maximum (hBest-like), as shipped since the original
              # implementation (#205).
              structure(charInfo / charMax, hMax = charMax)
@@ -356,16 +342,16 @@ ClusteringConcordance <- function(
            idx <- cbind(1, bestMatch, seq_along(bestMatch))
            return(weighted.mean(norm[idx], hBest[idx]))
           
-           if (isFALSE(normalize)) {
+           if (isFALSE(chanceCorrect)) {
            } else {
              one <- sum(hh["hChar", 1, ])
-             zero <- if (isTRUE(normalize)) {
+             zero <- if (isTRUE(chanceCorrect)) {
                sum(apply(hh["miRand", , ], 2, max))
              } else {
                randTreeMean
              }
              ret <- (sum(charInfo) - zero) / (one - zero)
-             if (is.numeric(normalize)) {
+             if (is.numeric(chanceCorrect)) {
                mcseInfo <- ((one - sum(charInfo)) / (one - zero) ^ 2) * mcseTree
                mcseInfo[mcseInfo < sqrt(.Machine$double.eps)] <- 0
                structure(ret, mcse = mcseInfo)
@@ -540,11 +526,11 @@ QALegend <- function(where = c(0.1, 0.3, 0.1, 0.3), n = 5, Col = QACol,
 #' @export
 ConcordanceTable <- function(tree, dataset, Col = QACol, largeClade = 0,
                              xlab = "Edge", ylab = "Character",
-                             normalize = TRUE, plot = TRUE,
+                             chanceCorrect = TRUE, plot = TRUE,
                              marginSize = 0L, paintSize = 0L,
                              palette = "default", ...) {
   cc <- ClusteringConcordance(tree, dataset, return = "all",
-                              normalize = normalize)
+                              chanceCorrect = chanceCorrect)
   nodes <- seq_len(dim(cc)[[2]])
   info <- .ConcSlice(cc, "hBest") * .ConcSlice(cc, "n")
   amount <- info / max(info, na.rm = TRUE)
@@ -743,9 +729,9 @@ MutualClusteringConcordance <- function(tree, dataset) {
 
 #' @rdname SiteConcordance
 #' @details
-#' `QuartetConcordance()` is the proportion of quartets (sets of four leaves)
-#' that are decisive for a split which are also concordant with it
-#' For example, a quartet with the characters `0 0 0 1` is not decisive, as
+#' `QuartetConcordance()` measures the agreement between each character and each
+#' split in the currency of quartets (sets of four leaves).
+#' A quartet with the characters `0 0 0 1` is not decisive, as
 #' all relationships between those leaves are equally parsimonious.
 #' But a quartet with characters `0 0 1 1` is decisive, and is concordant
 #' with any tree that groups the first two leaves together to the exclusion
@@ -755,12 +741,21 @@ MutualClusteringConcordance <- function(tree, dataset) {
 #' quartets that are decisive for a branch.
 #' Doing so circumvents the criticisms of \insertCite{Goloboff2024;textual}{TreeSearch}.
 #'
-#' By default, the reported value weights each site by the number of quartets
-#' it is decisive for.  This value can be interpreted as the proportion of
-#' all decisive quartets that are concordant with a split.
+#' The quartets that a split resolves are not logically independent: because
+#' \eqn{(ab, cd)} and \eqn{(ab, ce)} jointly entail \eqn{(ab, de)}
+#' \insertCite{Nelson1992}{TreeSearch}, a split of sizes \eqn{(k, t - k)}
+#' resolves \eqn{\binom{k}{2}\binom{t - k}{2}} quartets, of which only
+#' \eqn{(k - 1)(t - k - 1)} are non-redundant.
+#' `unit` selects which of the two is counted.  By default, agreement is
+#' measured against the non-redundant content, so that a character scores full
+#' marks only where it displays the split, rather than merely agreeing with a
+#' large combinatorial volume of its quartets.
+#'
+#' By default, the reported value weights each site by the quartet content it
+#' shares with the split.
 #' If `weight = FALSE`, the reported value is the mean of the concordance
-#' value for each site. 
-#' Consider a split associated with two sites:
+#' value for each site.
+#' Consider a split associated with two sites (counting in `unit = "quartet"`):
 #' one that is concordant with 25&percnt; of 96 decisive quartets, and
 #' a second that is concordant with 75&percnt; of 4 decisive quartets.
 #' If `weight = TRUE`, the split concordance will be 24 + 3 / 96 + 4 = 27&percnt;.
@@ -784,8 +779,20 @@ MutualClusteringConcordance <- function(tree, dataset) {
 #' Unlike the `"edge"` result, this vector is unnamed: characters retain no
 #' persistent identifier once collapsed to patterns.
 #'
+#' With `unit = "nrqs"`, per-character values are much smaller than 1 even
+#' where a character perfectly matches a split in a tree. Values rank characters
+#' on a given tree, but are not comparable between trees.
+#'
 #' @param weight Logical specifying whether to weight sites according to the
-#' number of quartets they are decisive for.
+#' quartet content that they share with each split.
+#' @param unit Character specifying the currency in which quartets are counted:
+#'   - `"nrqs"` (default): non-redundant quartet statements.
+#'     At a given split, only a character whose own split is identical scores
+#'     full marks; nested (compatible) characters receive partial support, and 
+#'     incompatible characters score lower still.
+#'   - `"quartet"`: each resolved quartet counts once,
+#'     analogous to the site concordance factor \insertCite{Minh2020}{TreeSearch}.
+#' @references \insertAllCited{}
 #' @importFrom ape keep.tip
 #' @importFrom cli cli_progress_bar cli_progress_update
 #' @importFrom utils combn
@@ -795,7 +802,9 @@ QuartetConcordance <- function(
   tree,
   dataset = NULL,
   weight = TRUE,
-  return = "edge"
+  return = "edge",
+  unit = c("nrqs", "quartet"),
+  chanceCorrect = TRUE
 ) {
   if (is.null(dataset)) {
     warning("Cannot calculate concordance without `dataset`.")
@@ -804,6 +813,15 @@ QuartetConcordance <- function(
   if (!inherits(dataset, "phyDat")) {
     stop("`dataset` must be a phyDat object.")
   }
+  unit <- match.arg(unit)
+  if (!isFALSE(chanceCorrect)) {
+    if (!isTRUE(chanceCorrect)) {
+      if (!is.numeric(chanceCorrect) || length(chanceCorrect) != 1L ||
+          is.na(chanceCorrect) || chanceCorrect < 1) {
+        stop("`chanceCorrect` must be FALSE, TRUE, or a positive integer.")
+      }
+    }
+  }
   tipLabels <- intersect(TipLabels(tree), names(dataset))
   if (!length(tipLabels)) {
     warning("No overlap between tree labels and dataset.")
@@ -840,10 +858,6 @@ QuartetConcordance <- function(
     dimnames = dimnames(characters)
   )
   
-  raw_counts <- quartet_concordance(logiSplits, charInt)
-
-  num <- raw_counts$concordant
-  den <- raw_counts$decisive
   options <- c("edge", "char")
   matched <- pmatch(tolower(trimws(return)), options, nomatch = NA_integer_)
   if (is.na(matched)) {
@@ -852,6 +866,34 @@ QuartetConcordance <- function(
   }
   return <- options[[matched]]
 
+  if (unit == "nrqs") {
+    # Return:
+    return(.NrqsConcordance(logiSplits, charInt, weight, return, splits,
+                            chanceCorrect))
+  }
+
+  raw_counts <- quartet_concordance(logiSplits, charInt)
+
+  num <- raw_counts$concordant
+  den <- raw_counts$decisive
+
+  # Chance correction (option): re-zero the observed concordant/decisive ratio
+  # against the ratio expected under the same fixed-marginal null used for the
+  # NRQS currency.  Only `conc` and `dec` vary under the null (the split sizes
+  # and state counts are fixed), so we need E[conc] and E[dec] per (split, char),
+  # computed exactly from the hypergeometric pmf (no floors here, unlike NRQS)
+  # or by Monte-Carlo tip-shuffle, then re-zero the pooled ratio.
+  doNorm <- !isFALSE(chanceCorrect)
+  if (doNorm) {
+    base <- if (isTRUE(chanceCorrect)) {
+      .QuartetExpect(charInt, logiSplits)
+    } else {
+      .QuartetMC(charInt, logiSplits, chanceCorrect)
+    }
+    eNum <- base[["concordant"]]
+    eDen <- base[["decisive"]]
+  }
+
   if (return == "edge") {
     if (isTRUE(weight)) {
       # Sum numerator and denominator across sites (columns), then divide
@@ -863,38 +905,360 @@ QuartetConcordance <- function(
         NA_real_,
         split_sums_num / split_sums_den
       )
+      if (doNorm) {
+        bDen <- rowSums(eDen)
+        base <- ifelse(bDen == 0, NA_real_, rowSums(eNum) / bDen)
+        ret <- .RezeroGuarded(ret, base)
+      }
     } else {
       # Mean of ratios per site
       # Avoid division by zero (0/0 -> NaN -> NA handled by na.rm)
       ratios <- num / den
       # Replace NaN/Inf with NA for rowMeans calculation
       ratios[!is.finite(ratios)] <- NA
-      ret <- rowMeans(ratios, na.rm = TRUE)
+      if (doNorm) {
+        bRatios <- eNum / eDen
+        bRatios[!is.finite(bRatios)] <- NA
+        # Average observed and baseline over the same (split, char) cells.
+        naMask <- is.na(ratios) | is.na(bRatios)
+        ratios[naMask] <- NA
+        bRatios[naMask] <- NA
+        ret <- .RezeroGuarded(rowMeans(ratios, na.rm = TRUE),
+                              rowMeans(bRatios, na.rm = TRUE))
+      } else {
+        ret <- rowMeans(ratios, na.rm = TRUE)
+      }
     }
 
     setNames(ret, names(splits))
   } else {
     # return = "char"
     if (isTRUE(weight)) {
-      vapply(
-        seq_len(dim(num)[[2]]),
-        function(i) {
-          weighted.mean(num[, i] / den[, i], den[, i])
-        },
-        double(1)
-      )
+      if (doNorm) {
+        obs <- vapply(seq_len(ncol(num)), function(i) {
+          d <- sum(den[, i]); if (d == 0) NA_real_ else sum(num[, i]) / d
+        }, double(1))
+        b <- vapply(seq_len(ncol(eNum)), function(i) {
+          d <- sum(eDen[, i]); if (d == 0) NA_real_ else sum(eNum[, i]) / d
+        }, double(1))
+        .RezeroGuarded(obs, b)
+      } else {
+        vapply(
+          seq_len(dim(num)[[2]]),
+          function(i) {
+            weighted.mean(num[, i] / den[, i], den[, i])
+          },
+          double(1)
+        )
+      }
     } else {
-      vapply(
-        seq_len(dim(num)[[2]]),
-        function(i) {
-          mean(num[den[, i] > 0, i] / den[den[, i] > 0, i])
-        },
-        double(1)
-      )
+      if (doNorm) {
+        sObs <- ifelse(den > 0, num / den, NA_real_)
+        sBase <- ifelse(eDen > 0, eNum / eDen, NA_real_)
+        naMask <- is.na(sObs) | is.na(sBase)
+        sObs[naMask] <- NA_real_
+        sBase[naMask] <- NA_real_
+        obs <- vapply(seq_len(ncol(sObs)), function(i) {
+          m <- mean(sObs[, i], na.rm = TRUE); if (is.nan(m)) NA_real_ else m
+        }, double(1))
+        b <- vapply(seq_len(ncol(sBase)), function(i) {
+          m <- mean(sBase[, i], na.rm = TRUE); if (is.nan(m)) NA_real_ else m
+        }, double(1))
+        .RezeroGuarded(obs, b)
+      } else {
+        vapply(
+          seq_len(dim(num)[[2]]),
+          function(i) {
+            mean(num[den[, i] > 0, i] / den[den[, i] > 0, i])
+          },
+          double(1)
+        )
+      }
     }
   }
 }
 
+# Expected raw concordant + decisive quartet counts under the fixed-marginal
+# null.  The exact expectation (per state-pair, an exact sum over the
+# trivariate-hypergeometric pmf) is computed in C++ (`quartet_expect` in
+# src/concordance_expect.cpp) for speed; `.QuartetExpect` is a thin R wrapper.
+# The `.QuartetMC` Monte-Carlo baseline reshuffles tokens and re-scores through
+# the C++ `quartet_concordance` kernel.  (The C++ path is validated bit-for-bit
+# against the earlier R implementation; see dev/benchmarks/frac-quart.)
+.QuartetExpect <- function(charInt, logiSplits) {
+  quartet_expect(logiSplits, charInt)
+}
+
+.QuartetMC <- function(charInt, logiSplits, nRelabel) {
+  nSplit <- ncol(logiSplits)
+  nChar <- ncol(charInt)
+  accConc <- accDec <- matrix(0, nSplit, nChar)
+  for (s in seq_len(nRelabel)) {
+    shuffled <- charInt
+    for (ci in seq_len(nChar)) {
+      col <- charInt[, ci]
+      scored <- !is.na(col)
+      shuffled[scored, ci] <- sample(col[scored])
+    }
+    kc <- quartet_concordance(logiSplits, shuffled)
+    accConc <- accConc + kc[["concordant"]]
+    accDec <- accDec + kc[["decisive"]]
+  }
+  list(concordant = accConc / nRelabel, decisive = accDec / nRelabel)
+}
+
+# Nelson-Ladiges fractional ("nrqs") currency for QuartetConcordance().
+#
+# The quartets a bipartition of sizes (k, t - k) resolves are the 4-cycles of
+# the complete bipartite graph K_{k, t-k}; the Nelson-Ladiges entailment
+# (ab,cd) + (ab,ce) -> (ab,de) is GF(2) cycle addition, so only the cyclomatic
+# number (k - 1)(t - k - 1) of them are non-redundant (the NRQS).
+#
+# A multi-state character is the disjoint union of its state-pairs: a quartet is
+# decisive only when two taxa share one state and two share another, so every
+# decisive quartet lives in the K_{n_i, n_j} block between two states.  Those
+# blocks are edge-disjoint and the entailment never crosses them (it forces the
+# shared "c, d, e" taxa into a single state), so NRQS add over state pairs:
+#   W_char = sum_{i<j} (n_i - 1)(n_j - 1),   A = sum_{i<j} A_ij
+# with a per state-pair weight of 4 / (n_i n_j).  (Verified by GF(2) rank of the
+# decisive / concordant quartet sets.)  Each state-pair is therefore an
+# independent binary sub-problem; a binary character is the single-pair case.
+#
+# Scoring uses coverage ("option b"): per pair, concordant NRQS over the
+# *reported* unit's own content (Wk for edges, Wc for characters), so that
+# A_ij <= min(Wc, Wk) keeps each ratio in [0, 1] and only an identical
+# character-split *pair* scores 1.
+#
+# That per-pair bound does NOT carry over to the value reported for
+# `return = "char"`.  The character score is the M-weighted mean of those
+# per-pair ratios across every split in the tree, and a character agrees
+# exactly with at most one split while merely being *compatible* with (silent
+# about) the rest.  Coverage scores silence as failure to cover, so even a
+# character identical to one of the tree's own splits averages well below 1:
+# 0.49 for a 22|26 split of a 48-leaf tree, 0.25 for a 2|46 split.  The
+# attainable maximum therefore depends on the tree, and character scores are
+# comparable to one another but not to an absolute ceiling of 1.
+#
+# There is no cheap repair.  Substituting the character's concordant NRQS with
+# respect to the whole tree would require counting a union of per-split
+# concordant sets, and those sets overlap heavily: summing A_ij across the 45
+# splits of the reference tree used in testing overshoots Wc by ~11x, so the
+# whole-tree count is not additive over splits and has no closed form here.
+# Taking max_e A_ij / Wc instead does attain 1, but best-match selection is
+# upward-biased in the absence of signal (see the `tree` return, which carries
+# the same warning) and E[max] is not available from the analytic per-split
+# expectations in `nrqs_expect()`, so `chanceCorrect = TRUE` could not support
+# it.  Callers who need a per-character score bounded at 1 should use
+# `unit = "quartet"`, whose char path reaches 1 for an identical character.
+#
+# Pairs are pooled by the shared information
+# M = min(Wc, Wk) (the hBest analogue), which is symmetric between character and
+# split so both `return`s pool by the same amount.  Ambiguous / inapplicable /
+# absent tokens drop out per character (treated as "?", as in the quartet path),
+# giving each character its own effective taxon count.
+.NrqsConcordance <- function(logiSplits, charInt, weight, return, splits,
+                             chanceCorrect) {
+  nSplit <- ncol(logiSplits)
+  nChar <- ncol(charInt)
+  pos <- function(z) {
+    z[z < 0] <- 0
+    z
+  }
+
+  # Observed contributions to the pools, summed over each character's
+  # state-pairs; the chance baseline (if requested) uses the same accumulators
+  # for the expected pools, so the re-zero divides like against like.
+  numEdge <- numChar <- denM <- matrix(0, nSplit, nChar)
+  wcTot <- numeric(nChar)                  # character NRQS content, split-free
+  doNorm <- !isFALSE(chanceCorrect)
+  if (doNorm) {
+    baseNumEdge <- baseNumChar <- baseDenM <- matrix(0, nSplit, nChar)
+    if (isTRUE(chanceCorrect)) {
+      # Exact hypergeometric expectation for all characters at once (C++;
+      # `nrqs_expect` in src/concordance_expect.cpp).  Uninformative characters
+      # yield zero pools, matching the observed loop's `wcTot > 0` guard.
+      be <- nrqs_expect(logiSplits, charInt)
+      baseNumEdge <- be[["numEdge"]]
+      baseNumChar <- be[["numChar"]]
+      baseDenM <- be[["denM"]]
+    }
+  }
+
+  for (ci in seq_len(nChar)) {
+    col <- charInt[, ci]
+    obs <- .CharNrqsContrib(col, logiSplits, pos)
+    numEdge[, ci] <- obs[["numEdge"]]
+    numChar[, ci] <- obs[["numChar"]]
+    denM[, ci] <- obs[["denM"]]
+    wcTot[ci] <- obs[["wcTot"]]
+    if (doNorm && !isTRUE(chanceCorrect) && obs[["wcTot"]] > 0) {
+      # Monte-Carlo shuffle
+      base <- .CharNrqsMC(col, logiSplits, pos, chanceCorrect)
+      baseNumEdge[, ci] <- base[["numEdge"]]
+      baseNumChar[, ci] <- base[["numChar"]]
+      baseDenM[, ci] <- base[["denM"]]
+    }
+  }
+
+  informative <- wcTot > 0
+
+  if (return == "edge") {
+    # edge: one value per split
+    if (isTRUE(weight)) {
+      denom <- rowSums(denM)
+      ret <- ifelse(denom == 0, NA_real_, rowSums(numEdge) / denom)
+      if (doNorm) {
+        bDen <- rowSums(baseDenM)
+        base <- ifelse(bDen == 0, NA_real_, rowSums(baseNumEdge) / bDen)
+        ret <- .RezeroGuarded(ret, base)
+      }
+    } else {
+      # Mean per-site quality over informative characters; uninformative
+      # characters carry no NRQS and are dropped, as in the quartet path.
+      sEdge <- ifelse(denM > 0, numEdge / denM, NA_real_)
+      if (doNorm) {
+        sBase <- ifelse(baseDenM > 0, baseNumEdge / baseDenM, NA_real_)
+        # Average observed and baseline over the SAME (split, char) cells: a
+        # degenerate multi-state pair can be dropped from one but not the other
+        # (observed wk = 0 yet E[m] > 0), which would otherwise re-zero mismatched
+        # populations.
+        naMask <- is.na(sEdge) | is.na(sBase)
+        sEdge[naMask] <- NA_real_
+        sBase[naMask] <- NA_real_
+        ret <- .RezeroGuarded(.RowMeanInformative(sEdge, informative, nSplit),
+                              .RowMeanInformative(sBase, informative, nSplit))
+      } else {
+        ret <- .RowMeanInformative(sEdge, informative, nSplit)
+      }
+    }
+    setNames(ret, names(splits))
+  } else {
+    # char: one value per character
+    if (isTRUE(weight)) {
+      denom <- colSums(denM)
+      ret <- ifelse(denom == 0, NA_real_, colSums(numChar) / denom)
+      if (doNorm) {
+        bDen <- colSums(baseDenM)
+        base <- ifelse(bDen == 0, NA_real_, colSums(baseNumChar) / bDen)
+        ret <- .RezeroGuarded(ret, base)
+      }
+      ret
+    } else {
+      sChar <- ifelse(denM > 0, numChar / denM, NA_real_)
+      if (doNorm) {
+        sBase <- ifelse(baseDenM > 0, baseNumChar / baseDenM, NA_real_)
+        # Match cells, as in the edge return above.
+        naMask <- is.na(sChar) | is.na(sBase)
+        sChar[naMask] <- NA_real_
+        sBase[naMask] <- NA_real_
+        ret <- .RezeroGuarded(.ColMeanInformative(sChar, informative, nChar),
+                              .ColMeanInformative(sBase, informative, nChar))
+      } else {
+        ret <- .ColMeanInformative(sChar, informative, nChar)
+      }
+      ret
+    }
+  }
+}
+
+# Observed per-character NRQS contributions, summed over the character's
+# state-pairs, returned as per-split vectors.  Shared by the observed pass and
+# the Monte-Carlo baseline so the two use byte-identical arithmetic.
+.CharNrqsContrib <- function(col, logiSplits, pos) {
+  nSplit <- ncol(logiSplits)
+  numEdge <- numChar <- denM <- numeric(nSplit)
+  wcTot <- 0
+  scored <- !is.na(col)
+  states <- sort(unique(col[scored]))
+  nStates <- length(states)
+  if (nStates >= 2L) {
+    for (a in seq_len(nStates - 1L)) {
+      for (b in seq(a + 1L, nStates)) {
+        inI <- scored & col == states[a]
+        inJ <- scored & col == states[b]
+        aI <- colSums(logiSplits & inI)      # state i, side A
+        bI <- colSums(!logiSplits & inI)     # state i, side B
+        aJ <- colSums(logiSplits & inJ)      # state j, side A
+        bJ <- colSums(!logiSplits & inJ)     # state j, side B
+        nI <- aI + bI                        # n_i (constant across splits)
+        nJ <- aJ + bJ                        # n_j
+        mA <- aI + aJ                        # taxa of this pair on side A
+        tP <- nI + nJ                        # taxa scored in this pair
+        # Concordant NRQS; the (x - 1)_+ floors stop self-agreement exceeding 1.
+        aij <- pos(aI - 1) * pos(bJ - 1) + pos(bI - 1) * pos(aJ - 1)
+        wc <- pos(nI - 1) * pos(nJ - 1)      # pair's character content
+        wk <- pos(mA - 1) * pos(tP - mA - 1) # pair's split content
+        m <- pmin(wc, wk)                    # shared information
+        denM <- denM + m
+        numEdge <- numEdge + ifelse(wk > 0, m * aij / wk, 0)
+        numChar <- numChar + ifelse(wc > 0, m * aij / wc, 0)
+        wcTot <- wcTot + wc[1]               # wc is constant across splits
+      }
+    }
+  }
+  list(numEdge = numEdge, numChar = numChar, denM = denM, wcTot = wcTot)
+}
+
+# The exact NRQS expectation E[m], E[m*A/wk], E[m*A/wc] is computed in C++
+# (`nrqs_expect` in src/concordance_expect.cpp) for all characters at once, and
+# consumed directly by `.NrqsConcordance`; there is no per-character R exact
+# helper.  The Monte-Carlo baseline below is the opt-in
+# `chanceCorrect = <int>` path.
+
+# Monte-Carlo baseline: average `.CharNrqsContrib` over `nRelabel` random
+# reassignments of the character's tokens across its scored leaves (the split,
+# the scored set and the state counts are held fixed).  Averaging the pools
+# (rather than per-split ratios) mirrors the exact ratio-of-expectations path.
+.CharNrqsMC <- function(col, logiSplits, pos, nRelabel) {
+  nSplit <- ncol(logiSplits)
+  scored <- !is.na(col)
+  tokens <- col[scored]
+  accEdge <- accChar <- accDen <- numeric(nSplit)
+  for (i in seq_len(nRelabel)) {
+    shuffled <- col
+    shuffled[scored] <- sample(tokens)
+    cc <- .CharNrqsContrib(shuffled, logiSplits, pos)
+    accEdge <- accEdge + cc[["numEdge"]]
+    accChar <- accChar + cc[["numChar"]]
+    accDen <- accDen + cc[["denM"]]
+  }
+  list(numEdge = accEdge / nRelabel,
+       numChar = accChar / nRelabel,
+       denM = accDen / nRelabel)
+}
+
+# Mean of per-split quality across informative characters, matching the
+# quartet path's treatment of uninformative characters (dropped).
+.RowMeanInformative <- function(mat, informative, nSplit) {
+  ret <- if (any(informative)) {
+    rowMeans(mat[, informative, drop = FALSE], na.rm = TRUE)
+  } else {
+    rep(NA_real_, nSplit)
+  }
+  ret[is.nan(ret)] <- NA_real_
+  ret
+}
+
+.ColMeanInformative <- function(mat, informative, nChar) {
+  vapply(seq_len(nChar), function(ci) {
+    if (informative[ci]) {
+      m <- mean(mat[, ci], na.rm = TRUE)
+      if (is.nan(m)) NA_real_ else m
+    } else {
+      NA_real_
+    }
+  }, double(1))
+}
+
+# Re-zero to a chance baseline, guarding the degenerate `zero -> 1` case
+# (returns NA rather than dividing by ~0).  Values are returned UNCLAMPED, as
+# in `ClusteringConcordance`; clamping to [-1, 1] is deferred to plotting.
+.RezeroGuarded <- function(value, zero) {
+  denom <- 1 - zero
+  ifelse(is.na(value) | is.na(zero) | denom < sqrt(.Machine[["double.eps"]]),
+         NA_real_, (value - zero) / denom)
+}
+
 .ExpectedMICache <- new.env(hash = TRUE, parent = emptyenv())
 # Bound on cache size: a session computing concordance for many differently
 # -sized trees/characters would otherwise grow this cache without limit.
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 7e94f7957..348ed087b 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -20,6 +20,14 @@ MaddisonSlatkin_clear_cache <- function() {
     invisible(.Call(`_TreeSearch_MaddisonSlatkin_clear_cache`))
 }
 
+quartet_expect <- function(splits, characters) {
+    .Call(`_TreeSearch_quartet_expect`, splits, characters)
+}
+
+nrqs_expect <- function(splits, characters) {
+    .Call(`_TreeSearch_nrqs_expect`, splits, characters)
+}
+
 #' Expected mutual information between two partitions
 #'
 #' Computes the mutual information expected purely by chance between two
@@ -27,7 +35,7 @@ MaddisonSlatkin_clear_cache <- function() {
 #' block sizes (marginals) of each partition are fixed but the items are
 #' associated at random.  Subtracting this baseline from an observed mutual
 #' information yields a chance-corrected ("adjusted") mutual information, as
-#' applied by [`SiteConcordance`]`(normalize = TRUE)`; it is most material for
+#' applied by [`SiteConcordance`]`(chanceCorrect = TRUE)`; it is most material for
 #' small or unbalanced partitions, where raw mutual information is appreciably
 #' inflated by chance agreement.
 #'
diff --git a/dev/benchmarks/frac-quart/chance_baseline.R b/dev/benchmarks/frac-quart/chance_baseline.R
new file mode 100644
index 000000000..38187da82
--- /dev/null
+++ b/dev/benchmarks/frac-quart/chance_baseline.R
@@ -0,0 +1,86 @@
+# Validate the exact hypergeometric chance baseline used by
+# QuartetConcordance(unit = "nrqs", chanceCorrect = TRUE) against a Monte-Carlo
+# tip-shuffle oracle, at the level of a single (split, character, state-pair).
+#
+# Under the fixed-marginal null the character's tokens are reassigned across the
+# t scored leaves at random, holding the pair's state counts (n_i, n_j) and the
+# split's scored side-A size M fixed.  The cell vector (p,q,r,s) is then
+# trivariate-hypergeometric.  We need the EXPECTED per-pair pool contributions
+#   m       = min(w_c, w_k)          (shared-information pooling weight)
+#   m*A/w_k (edge-numerator contribution)
+#   m*A/w_c (char-numerator contribution)
+# because w_k = (mA-1)+(tP-mA-1)+ -- and hence m -- are THEMSELVES random for
+# multistate pairs (mA = p + r varies).  Binary characters are the special case
+# mA == M, tP == t, so w_k is fixed and only A varies (the clean
+# ClusteringConcordance case).  This mirrors `.ExpectedTrit()` in R/Concordance.R.
+
+pos <- function(z) { z[z < 0] <- 0; z }
+
+# Observed per-pair triple for one realised cell vector.
+pairTriple <- function(p, q, r, s) {
+  nI <- p + q; nJ <- r + s
+  mA <- p + r; tP <- nI + nJ
+  A  <- pos(p - 1) * pos(s - 1) + pos(q - 1) * pos(r - 1)
+  wc <- pos(nI - 1) * pos(nJ - 1)
+  wk <- pos(mA - 1) * pos(tP - mA - 1)
+  m  <- min(wc, wk)
+  c(m    = m,
+    mAwk = if (wk > 0) m * A / wk else 0,
+    mAwc = if (wc > 0) m * A / wc else 0)
+}
+
+# Exact expectation over the trivariate hypergeometric (the R kernel's method).
+expectedTriple <- function(nI, nJ, M, t) {
+  nOther <- t - nI - nJ
+  stopifnot(nOther >= 0, M >= 0, M <= t)
+  logC <- function(n, k) if (k < 0 || k > n) -Inf else lchoose(n, k)
+  denom <- lchoose(t, M)
+  acc <- c(m = 0, mAwk = 0, mAwc = 0)
+  for (p in 0:min(nI, M)) {
+    for (r in 0:min(nJ, M - p)) {
+      oA <- M - p - r
+      if (oA < 0 || oA > nOther) next
+      prob <- exp(logC(nI, p) + logC(nJ, r) + logC(nOther, oA) - denom)
+      if (prob <= 0) next
+      acc <- acc + prob * pairTriple(p, nI - p, r, nJ - r)
+    }
+  }
+  acc
+}
+
+# Monte-Carlo oracle: shuffle state labels across t leaves, M fixed on side A.
+mcTriple <- function(nI, nJ, M, t, nShuf = 2e5) {
+  sideA <- c(rep(TRUE, M), rep(FALSE, t - M))
+  acc <- c(m = 0, mAwk = 0, mAwc = 0)
+  for (i in seq_len(nShuf)) {
+    lab <- sample(rep(1:3, c(nI, nJ, t - nI - nJ)))   # 1 = i, 2 = j, 3 = other
+    p <- sum(lab == 1 & sideA); q <- sum(lab == 1 & !sideA)
+    r <- sum(lab == 2 & sideA); s <- sum(lab == 2 & !sideA)
+    acc <- acc + pairTriple(p, q, r, s)
+  }
+  acc / nShuf
+}
+
+set.seed(1)
+cases <- list(
+  c(nI = 4, nJ = 4, M = 4, t = 8),   # binary-like (nOther = 0): w_k fixed
+  c(nI = 3, nJ = 5, M = 4, t = 8),
+  c(nI = 3, nJ = 3, M = 4, t = 9),   # multistate (nOther > 0): w_k random
+  c(nI = 4, nJ = 2, M = 5, t = 10),
+  c(nI = 5, nJ = 4, M = 6, t = 12),
+  c(nI = 2, nJ = 2, M = 3, t = 7),
+  c(nI = 6, nJ = 3, M = 5, t = 15)
+)
+
+worst <- 0
+for (cs in cases) {
+  ex <- expectedTriple(cs["nI"], cs["nJ"], cs["M"], cs["t"])
+  mc <- mcTriple(cs["nI"], cs["nJ"], cs["M"], cs["t"])
+  worst <- max(worst, max(abs(ex - mc)))
+  cat(sprintf("nI=%d nJ=%d M=%2d t=%2d  exact=(%.4f,%.4f,%.4f)  mc=(%.4f,%.4f,%.4f)\n",
+              cs["nI"], cs["nJ"], cs["M"], cs["t"],
+              ex[1], ex[2], ex[3], mc[1], mc[2], mc[3]))
+}
+cat(sprintf("\nWorst exact-vs-MC discrepancy: %.5f\n", worst))
+stopifnot(worst < 0.05)
+cat("PASS: exact hypergeometric baseline matches the MC tip-shuffle oracle.\n")
diff --git a/dev/benchmarks/frac-quart/cpp_expect_parity.R b/dev/benchmarks/frac-quart/cpp_expect_parity.R
new file mode 100644
index 000000000..10972edd3
--- /dev/null
+++ b/dev/benchmarks/frac-quart/cpp_expect_parity.R
@@ -0,0 +1,108 @@
+# Parity guard for the C++ chance-baseline kernels `quartet_expect` and
+# `nrqs_expect` (src/concordance_expect.cpp), which compute the exact expected
+# concordant/decisive quartet counts (raw) and the expected trit pools
+# (E[m], E[m*A/wk], E[m*A/wc]) under the fixed-marginal (hypergeometric) null.
+#
+# The kernels replaced an earlier pure-R implementation; this script keeps a
+# self-contained R reference (the same trivariate-hypergeometric double sum) and
+# checks the C++ output matches it to machine precision across random binary,
+# multistate and missing-data matrices.  Run after any change to the kernels.
+
+suppressMessages(pkgload::load_all(
+  "C:/Users/pjjg18/GitHub/worktrees/TreeSearch/frac-quart-cpp",
+  quiet = TRUE, recompile = TRUE))
+suppressMessages(library(TreeTools))
+
+# Build logiSplits + charInt exactly as QuartetConcordance() does internally.
+build_inputs <- function(tree, dat) {
+  tl <- intersect(TipLabels(tree), names(dat))
+  dat <- dat[tl, drop = FALSE]
+  splits <- as.Splits(tree, dat)
+  logiSplits <- vapply(seq_along(splits), function(i) as.logical(splits[[i]]),
+                       logical(NTip(dat)))
+  contrast <- attr(dat, "contrast"); lv <- attr(dat, "allLevels")
+  isInapp <- lv == "-"
+  isAmbig <- rowSums(contrast[, colnames(contrast) != "-", drop = FALSE]) > 1
+  isGrouping <- !isAmbig & !isInapp
+  gCols <- apply(contrast[isGrouping, , drop = FALSE] > 0, 1, which)
+  l2i <- rep(NA_integer_, length(lv)); l2i[isGrouping] <- as.integer(gCols)
+  ch <- PhyDatToMatrix(dat)
+  charInt <- array(l2i[match(ch, lv)], dim = dim(ch), dimnames = dimnames(ch))
+  list(logiSplits = logiSplits, charInt = charInt)
+}
+
+# ---- self-contained R reference (trivariate-hypergeometric double sum) ----
+pos1 <- function(z) if (z < 0) 0 else z
+choose2 <- function(z) if (z < 2) 0 else z * (z - 1) / 2
+
+ref_pair <- function(nI, nJ, M, t) {           # returns c(Econc, Edec, Em, Ewk, Ewc)
+  nOther <- t - nI - nJ; tP <- nI + nJ
+  wc <- pos1(nI - 1) * pos1(nJ - 1)
+  logDen <- lchoose(t, M)
+  acc <- c(0, 0, 0, 0, 0)
+  for (p in 0:min(nI, M)) for (r in 0:min(nJ, M - p)) {
+    oA <- M - p - r
+    if (oA < 0 || oA > nOther) next
+    prob <- exp(lchoose(nI, p) + lchoose(nJ, r) + lchoose(nOther, oA) - logDen)
+    if (prob <= 0) next
+    q <- nI - p; s <- nJ - r; mA <- p + r
+    conc <- choose2(p) * choose2(s) + choose2(q) * choose2(r)
+    dec <- conc + p * q * r * s
+    A <- pos1(p - 1) * pos1(s - 1) + pos1(q - 1) * pos1(r - 1)
+    wk <- pos1(mA - 1) * pos1(tP - mA - 1); m <- min(wc, wk)
+    acc <- acc + prob * c(conc, dec, m, if (wk > 0) m * A / wk else 0,
+                          if (wc > 0) m * A / wc else 0)
+  }
+  acc
+}
+
+ref_expect <- function(charInt, logiSplits) {
+  nSplit <- ncol(logiSplits); nChar <- ncol(charInt)
+  out <- list(eConc = matrix(0, nSplit, nChar), eDec = matrix(0, nSplit, nChar),
+              eM = matrix(0, nSplit, nChar), eWk = matrix(0, nSplit, nChar),
+              eWc = matrix(0, nSplit, nChar))
+  for (ci in seq_len(nChar)) {
+    col <- charInt[, ci]; scored <- !is.na(col)
+    states <- sort(unique(col[scored])); nStates <- length(states)
+    if (nStates < 2L) next
+    tc <- sum(scored); mSideA <- colSums(logiSplits & scored)
+    uM <- unique(mSideA); idx <- match(mSideA, uM)
+    cnt <- tabulate(match(col[scored], states), nStates)
+    for (a in seq_len(nStates - 1L)) for (b in seq(a + 1L, nStates)) {
+      e <- vapply(uM, function(M) ref_pair(cnt[a], cnt[b], M, tc), double(5))
+      out$eConc[, ci] <- out$eConc[, ci] + e[1, idx]
+      out$eDec[, ci]  <- out$eDec[, ci]  + e[2, idx]
+      out$eM[, ci]    <- out$eM[, ci]    + e[3, idx]
+      out$eWk[, ci]   <- out$eWk[, ci]   + e[4, idx]
+      out$eWc[, ci]   <- out$eWc[, ci]   + e[5, idx]
+    }
+  }
+  out
+}
+
+set.seed(7)
+worst <- 0
+for (rep in 1:10) {
+  nTip <- sample(9:44, 1)
+  tree <- RandomTree(nTip, root = FALSE)
+  m <- replicate(sample(5:30, 1), {
+    ns <- sample(2:4, 1)
+    x <- as.character(sample(0:(ns - 1), nTip, replace = TRUE))
+    if (runif(1) < 0.3) x[sample(nTip, sample(1:3, 1))] <- "?"
+    x
+  })
+  rownames(m) <- paste0("t", seq_len(nTip))
+  io <- build_inputs(tree, MatrixToPhyDat(m))
+  ls <- io$logiSplits; ci <- io$charInt
+  R <- ref_expect(ci, ls)
+  cq <- quartet_expect(ls, ci); ct <- nrqs_expect(ls, ci)
+  d <- max(abs(cq$concordant - R$eConc), abs(cq$decisive - R$eDec),
+           abs(ct$numEdge - R$eWk), abs(ct$numChar - R$eWc),
+           abs(ct$denM - R$eM))
+  worst <- max(worst, d)
+  cat(sprintf("rep %2d: nTip=%2d nChar=%2d  max|C++ - R| = %.2e\n",
+              rep, nTip, ncol(ci), d))
+}
+cat(sprintf("\nWorst |C++ - R| over all kernels: %.3e\n", worst))
+stopifnot(worst < 1e-9)
+cat("PASS: C++ expectation kernels match the R reference to machine precision.\n")
diff --git a/dev/benchmarks/frac-quart/multistate_trit.R b/dev/benchmarks/frac-quart/multistate_trit.R
new file mode 100644
index 000000000..d11e38e1d
--- /dev/null
+++ b/dev/benchmarks/frac-quart/multistate_trit.R
@@ -0,0 +1,88 @@
+# Verify: for a multistate character, GF(2) rank of its decisive quartets
+# equals Sum_{i<j}(n_i-1)(n_j-1), and rank of concordant quartets (vs a split)
+# equals Sum_{i<j} A_ij.  Blocks are edge-disjoint => trits add across pairs.
+
+gf2rank <- function(M) {                      # rank over GF(2)
+  if (nrow(M) == 0 || ncol(M) == 0) return(0L)
+  M <- M %% 2L; r <- 0L; nc <- ncol(M)
+  for (col in seq_len(nc)) {
+    piv <- which(M[, col] == 1L)
+    piv <- piv[piv > r]
+    if (!length(piv)) next
+    p <- piv[1]; r <- r + 1L
+    M[c(r, p), ] <- M[c(p, r), ]
+    others <- which(M[, col] == 1L); others <- others[others != r]
+    if (length(others)) M[others, ] <- sweep(M[others, , drop = FALSE], 2, M[r, ], `+`) %% 2L
+    if (r == nrow(M)) break
+  }
+  r
+}
+
+# edge id for unordered pair (a,b) among t taxa
+eid <- function(a, b, t) { if (a > b) { tmp <- a; a <- b; b <- tmp }; (a - 1) * t + b }
+
+# All decisive quartets of `char` (state vector), each as its char-induced
+# 4-cycle edge vector over C(t,2)... encoded on t*t (sparse-safe) columns.
+quartet_rows <- function(char, split = NULL, concordant_only = FALSE) {
+  t <- length(char)
+  qs <- combn(t, 4)
+  rows <- list()
+  for (k in seq_len(ncol(qs))) {
+    q <- qs[, k]; st <- char[q]
+    tb <- table(st)
+    if (length(tb) == 2 && all(tb == 2)) {          # "xx yy": decisive
+      states <- as.integer(names(tb))
+      g1 <- q[st == states[1]]; g2 <- q[st == states[2]]   # the two same-state pairs
+      if (concordant_only) {
+        s1 <- split[g1]; s2 <- split[g2]
+        # concordant iff each same-state pair lies wholly on one split side,
+        # and the two pairs are on opposite sides
+        ok <- length(unique(s1)) == 1 && length(unique(s2)) == 1 && s1[1] != s2[1]
+        if (!ok) next
+      }
+      v <- integer(t * t)
+      for (a in g1) for (b in g2) v[eid(a, b, t)] <- 1L    # cross edges = 4-cycle
+      rows[[length(rows) + 1]] <- v
+    }
+  }
+  if (!length(rows)) return(matrix(0L, 0, t * t))
+  do.call(rbind, rows)
+}
+
+Aij <- function(p, q, r, s) max(p - 1, 0) * max(s - 1, 0) + max(q - 1, 0) * max(r - 1, 0)
+
+set.seed(1)
+cat(sprintf("%-28s %6s %6s %6s\n", "case", "rankD", "SumW", "match"))
+for (trial in 1:6) {
+  t <- sample(6:9, 1)
+  r <- sample(2:4, 1)
+  char <- sample(seq_len(r), t, replace = TRUE)
+  while (length(unique(char)) < 2 || any(tabulate(char) == 1 & tabulate(char) > 0 & FALSE)) char <- sample(seq_len(r), t, replace = TRUE)
+  ns <- as.integer(tabulate(char)); ns <- ns[ns > 0]
+  sumW <- sum(outer(seq_along(ns), seq_along(ns), Vectorize(function(i, j)
+    if (i < j) (ns[i] - 1) * (ns[j] - 1) else 0)))
+  rD <- gf2rank(quartet_rows(char))
+  cat(sprintf("t=%d states=%-14s %6d %6d %6s\n",
+              t, paste(ns, collapse = ","), rD, sumW, rD == sumW))
+}
+
+cat("\n-- concordant vs split --\n")
+cat(sprintf("%-34s %6s %6s %6s\n", "case", "rankC", "SumAij", "match"))
+for (trial in 1:6) {
+  t <- sample(6:9, 1); r <- sample(2:4, 1)
+  char <- sample(seq_len(r), t, replace = TRUE)
+  while (length(unique(char)) < 2) char <- sample(seq_len(r), t, replace = TRUE)
+  split <- sample(c(TRUE, FALSE), t, replace = TRUE)
+  states <- sort(unique(char))
+  sumA <- 0
+  for (ii in seq_along(states)) for (jj in seq_along(states)) if (ii < jj) {
+    i <- states[ii]; j <- states[jj]
+    p <- sum(char == i & split); qc <- sum(char == i & !split)
+    rr <- sum(char == j & split); ss <- sum(char == j & !split)
+    sumA <- sumA + Aij(p, qc, rr, ss)
+  }
+  rC <- gf2rank(quartet_rows(char, split, concordant_only = TRUE))
+  cat(sprintf("t=%d states=%-12s split=%-2d/%-2d %6d %6d %6s\n",
+              t, paste(as.integer(tabulate(char))[tabulate(char) > 0], collapse = ","),
+              sum(split), sum(!split), rC, sumA, rC == sumA))
+}
diff --git a/dev/benchmarks/frac-quart/validate_trit.R b/dev/benchmarks/frac-quart/validate_trit.R
new file mode 100644
index 000000000..b29e053c6
--- /dev/null
+++ b/dev/benchmarks/frac-quart/validate_trit.R
@@ -0,0 +1,229 @@
+# Validation of the unit = "nrqs" path in QuartetConcordance().
+#
+# Ordering follows the advisor's priority:
+#   1. Cell extraction cross-checked against the verified C++ kernel (conc/dec).
+#   2. A_ij <= min(Wc, Wk) asserted on every real (split, char, pair).
+#   3. Spec ladder (t = 8) reproduced from the per-pair formula.
+#   4. An INDEPENDENT re-implementation of the measure matches the package on
+#      binary (congreveLamsdell), multistate, and missing-data inputs.
+#   5. `return` parsing quirk mirrored between quartet and trit paths.
+#   6. congreveLamsdell quartet -> trit edge-score movement + regime census.
+#
+# Run: Rscript dev/benchmarks/frac-quart/validate_trit.R
+
+suppressMessages(pkgload::load_all(".", quiet = TRUE))
+library("TreeTools", quietly = TRUE)
+
+ok <- function(label, cond) {
+  cat(sprintf("[%s] %s\n", if (isTRUE(cond)) "PASS" else "FAIL", label))
+  if (!isTRUE(cond)) stop("FAILED: ", label)
+}
+
+# ---- Rebuild the (logiSplits, charInt) the function feeds to its kernels ----
+build_inputs <- function(tree, dataset) {
+  tipLabels <- intersect(TipLabels(tree), names(dataset))
+  dataset <- dataset[tipLabels, drop = FALSE]
+  splits <- as.Splits(tree, dataset)
+  logiSplits <- vapply(seq_along(splits), function(i) as.logical(splits[[i]]),
+                       logical(NTip(dataset)))
+  contrast <- attr(dataset, "contrast")
+  charLevels <- attr(dataset, "allLevels")
+  isInapp <- charLevels == "-"
+  isAmbig <- rowSums(contrast[, colnames(contrast) != "-", drop = FALSE]) > 1
+  isGrouping <- !isAmbig & !isInapp
+  groupingCols <- apply(contrast[isGrouping, , drop = FALSE] > 0, 1, which)
+  levelToInt <- rep(NA_integer_, length(charLevels))
+  levelToInt[isGrouping] <- as.integer(groupingCols)
+  characters <- PhyDatToMatrix(dataset)
+  charInt <- array(levelToInt[match(characters, charLevels)],
+                   dim = dim(characters), dimnames = dimnames(characters))
+  list(logiSplits = logiSplits, charInt = charInt, splits = splits)
+}
+
+ch2 <- function(x) x * (x - 1) / 2
+posm <- function(z) { z[z < 0] <- 0; z }
+
+# Per (split, char) cells summed over state-pairs, returned as component lists
+# so the same crosstab feeds both the kernel cross-check and the trit measure.
+pair_cells <- function(logiSplits, charInt) {
+  nSplit <- ncol(logiSplits); nChar <- ncol(charInt)
+  out <- vector("list", nChar)
+  for (ci in seq_len(nChar)) {
+    col <- charInt[, ci]; scored <- !is.na(col)
+    states <- sort(unique(col[scored]))
+    prs <- list()
+    if (length(states) >= 2) {
+      for (a in seq_len(length(states) - 1L)) for (b in seq(a + 1L, length(states))) {
+        inI <- scored & col == states[a]; inJ <- scored & col == states[b]
+        prs[[length(prs) + 1L]] <- list(
+          aI = colSums(logiSplits & inI), bI = colSums(!logiSplits & inI),
+          aJ = colSums(logiSplits & inJ), bJ = colSums(!logiSplits & inJ))
+      }
+    }
+    out[[ci]] <- prs
+  }
+  out
+}
+
+# ---------- 1. Cell extraction vs the verified kernel ----------
+{
+  data("congreveLamsdellMatrices", package = "TreeSearch")
+  dataset <- congreveLamsdellMatrices[[1]]
+  tree <- TreeSearch::referenceTree
+  inp <- build_inputs(tree, dataset)
+  cells <- pair_cells(inp$logiSplits, inp$charInt)
+  nSplit <- ncol(inp$logiSplits); nChar <- ncol(inp$charInt)
+  reconC <- reconD <- matrix(0, nSplit, nChar)
+  for (ci in seq_len(nChar)) for (p in cells[[ci]]) {
+    reconC[, ci] <- reconC[, ci] + ch2(p$aI) * ch2(p$bJ) + ch2(p$bI) * ch2(p$aJ)
+    reconD[, ci] <- reconD[, ci] + p$aI * p$bI * p$aJ * p$bJ
+  }
+  reconD <- reconD + reconC
+  kern <- TreeSearch:::quartet_concordance(inp$logiSplits, inp$charInt)
+  ok("cell -> concordant matches kernel", all(abs(reconC - kern$concordant) < 1e-9))
+  ok("cell -> decisive   matches kernel", all(abs(reconD - kern$decisive) < 1e-9))
+
+  # ---------- 2. A_ij <= min(Wc, Wk) on every pair ----------
+  worst <- -Inf
+  for (ci in seq_len(nChar)) for (p in cells[[ci]]) {
+    aij <- posm(p$aI - 1) * posm(p$bJ - 1) + posm(p$bI - 1) * posm(p$aJ - 1)
+    nI <- p$aI + p$bI; nJ <- p$aJ + p$bJ; mA <- p$aI + p$aJ; tP <- nI + nJ
+    wc <- posm(nI - 1) * posm(nJ - 1); wk <- posm(mA - 1) * posm(tP - mA - 1)
+    worst <- max(worst, max(aij - pmin(wc, wk)))
+  }
+  ok("A_ij <= min(Wc, Wk) everywhere", worst <= 1e-9)
+}
+
+# ---------- 3. Spec ladder (t = 8) from the per-pair formula ----------
+pairQ <- function(p, q, r, s) {           # edge quality A / Wk for one 2x2 table
+  A <- posm(p - 1) * posm(s - 1) + posm(q - 1) * posm(r - 1)
+  m <- p + r; t <- p + q + r + s
+  Wk <- posm(m - 1) * posm(t - m - 1)
+  A / Wk
+}
+ladder <- rbind(
+  identical    = c(4, 0, 0, 4),
+  nested       = c(3, 0, 2, 3),
+  mildCrossing = c(3, 1, 1, 3),
+  maxCrossing  = c(2, 2, 2, 2))
+Qlad <- apply(ladder, 1, function(x) pairQ(x[1], x[2], x[3], x[4]))
+expected <- c(identical = 1, nested = 0.5, mildCrossing = 4 / 9, maxCrossing = 2 / 9)
+ok("ladder reproduces spec table", all(abs(Qlad - expected) < 1e-9))
+ok("only identical scores 1", sum(abs(Qlad - 1) < 1e-9) == 1L)
+ok("crossing < nested", Qlad["mildCrossing"] < Qlad["nested"] &&
+     Qlad["maxCrossing"] < Qlad["mildCrossing"])
+
+# ---------- 4. Independent re-implementation vs the package ----------
+# Deliberately different code shape: per-pair matrices, explicit pooling.
+ref_trit <- function(tree, dataset, weight = TRUE, return = "edge") {
+  inp <- build_inputs(tree, dataset)
+  cells <- pair_cells(inp$logiSplits, inp$charInt)
+  nSplit <- ncol(inp$logiSplits); nChar <- ncol(inp$charInt)
+  numE <- numC <- den <- matrix(0, nSplit, nChar); wcTot <- numeric(nChar)
+  for (ci in seq_len(nChar)) for (p in cells[[ci]]) {
+    aij <- posm(p$aI - 1) * posm(p$bJ - 1) + posm(p$bI - 1) * posm(p$aJ - 1)
+    nI <- p$aI + p$bI; nJ <- p$aJ + p$bJ; mA <- p$aI + p$aJ; tP <- nI + nJ
+    wc <- posm(nI - 1) * posm(nJ - 1); wk <- posm(mA - 1) * posm(tP - mA - 1)
+    m <- pmin(wc, wk)
+    den[, ci] <- den[, ci] + m
+    numE[, ci] <- numE[, ci] + ifelse(wk > 0, m * aij / wk, 0)
+    numC[, ci] <- numC[, ci] + ifelse(wc > 0, m * aij / wc, 0)
+    wcTot[ci] <- wcTot[ci] + wc[1]
+  }
+  info <- wcTot > 0
+  ret <- pmatch(tolower(trimws(return)), c("character", "site", "default"),
+                nomatch = 3L)
+  if (ret == 3L) {                                  # edge
+    if (isTRUE(weight)) {
+      d <- rowSums(den); v <- ifelse(d == 0, NA_real_, rowSums(numE) / d)
+    } else {
+      sE <- ifelse(den > 0, numE / den, NA_real_)
+      v <- if (any(info)) rowMeans(sE[, info, drop = FALSE], na.rm = TRUE)
+           else rep(NA_real_, nSplit)
+      v[is.nan(v)] <- NA_real_
+    }
+    setNames(v, names(inp$splits))
+  } else {                                          # char
+    if (isTRUE(weight)) {
+      d <- colSums(den); ifelse(d == 0, NA_real_, colSums(numC) / d)
+    } else {
+      sC <- ifelse(den > 0, numC / den, NA_real_)
+      vapply(seq_len(nChar), function(ci) if (info[ci]) {
+        mm <- mean(sC[, ci], na.rm = TRUE); if (is.nan(mm)) NA_real_ else mm
+      } else NA_real_, double(1))
+    }
+  }
+}
+
+same <- function(a, b) all((is.na(a) & is.na(b)) | (abs(a - b) < 1e-9), na.rm = FALSE)
+
+data("congreveLamsdellMatrices", package = "TreeSearch")
+binDat <- congreveLamsdellMatrices[[1]]
+refTree <- TreeSearch::referenceTree
+for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) {
+  pkg <- QuartetConcordance(refTree, binDat, weight = w, return = r, unit = "nrqs", chanceCorrect = FALSE)
+  rf  <- ref_trit(refTree, binDat, weight = w, return = r)
+  ok(sprintf("binary  pkg==ref  weight=%s return=%s", w, r), same(pkg, rf))
+}
+
+# Multistate toy (3- and 4-state characters, complete)
+msMat <- matrix(c(0, 0, 1, 1, 2, 2, 2, 0,     # 3 states
+                  0, 0, 0, 1, 1, 2, 3, 3,     # 4 states
+                  0, 0, 1, 1, 2, 2, 0, 1), 8, # 3 states
+                dimnames = list(paste0("t", 1:8), NULL))
+msDat <- MatrixToPhyDat(msMat)
+msTree <- BalancedTree(8)
+for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) {
+  pkg <- QuartetConcordance(msTree, msDat, weight = w, return = r, unit = "nrqs", chanceCorrect = FALSE)
+  rf  <- ref_trit(msTree, msDat, weight = w, return = r)
+  ok(sprintf("multi   pkg==ref  weight=%s return=%s", w, r), same(pkg, rf))
+}
+
+# Missing / ambiguous / inapplicable toy
+naMat <- matrix(c(0, 0, 1, 1, "?", "?", 0, 1,
+                  0, 1, "-", 1, 0, "(01)", 1, 0,
+                  0, 0, 0, 1, 1, 2, "?", 2), 8,
+                dimnames = list(paste0("t", 1:8), NULL))
+naDat <- MatrixToPhyDat(naMat)
+for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) {
+  pkg <- QuartetConcordance(msTree, naDat, weight = w, return = r, unit = "nrqs", chanceCorrect = FALSE)
+  rf  <- ref_trit(msTree, naDat, weight = w, return = r)
+  ok(sprintf("missing pkg==ref  weight=%s return=%s", w, r), same(pkg, rf))
+}
+
+# ---------- 5. `return` parsing mirrors the quartet path ----------
+e1 <- QuartetConcordance(refTree, binDat, return = "edge",      unit = "nrqs", chanceCorrect = FALSE)
+e2 <- QuartetConcordance(refTree, binDat, return = "default",   unit = "nrqs", chanceCorrect = FALSE)
+c1 <- QuartetConcordance(refTree, binDat, return = "char",      unit = "nrqs", chanceCorrect = FALSE)
+c2 <- QuartetConcordance(refTree, binDat, return = "character", unit = "nrqs", chanceCorrect = FALSE)
+c3 <- QuartetConcordance(refTree, binDat, return = "site",      unit = "nrqs", chanceCorrect = FALSE)
+ok("return edge == default", same(e1, e2))
+ok("return char == character == site", same(c1, c2) && same(c1, c3))
+ok("edge length == n splits, char length == n chars",
+   length(e1) == length(inp <- as.Splits(refTree)) &&
+     length(c1) == sum(attr(binDat, "weight")))
+
+# ---------- 6. congreveLamsdell quartet -> trit movement + census ----------
+qEdge <- QuartetConcordance(refTree, binDat, unit = "quartet", chanceCorrect = FALSE)
+tEdge <- QuartetConcordance(refTree, binDat, unit = "nrqs", chanceCorrect = FALSE)
+cat("\nEdge concordance, quartet vs trit (congreveLamsdell[[1]]):\n")
+print(round(rbind(quartet = qEdge, trit = tEdge, delta = tEdge - qEdge), 3))
+cat(sprintf("\nmean quartet = %.3f  mean trit = %.3f  (trit is stricter)\n",
+            mean(qEdge, na.rm = TRUE), mean(tEdge, na.rm = TRUE)))
+
+# Regime census per edge: identical / nested / crossing state-pairs
+inpc <- build_inputs(refTree, binDat)
+cells <- pair_cells(inpc$logiSplits, inpc$charInt)
+regime <- matrix(0L, ncol(inpc$logiSplits), 3,
+                 dimnames = list(names(inpc$splits),
+                                 c("identical", "nested", "crossing")))
+for (ci in seq_along(cells)) for (p in cells[[ci]]) {
+  empt <- (p$aI == 0) + (p$bI == 0) + (p$aJ == 0) + (p$bJ == 0)
+  regime[, "identical"] <- regime[, "identical"] + (empt == 2)
+  regime[, "nested"]    <- regime[, "nested"]    + (empt == 1)
+  regime[, "crossing"]  <- regime[, "crossing"]  + (empt == 0)
+}
+cat("\nState-pair regime census per edge:\n")
+print(regime)
+
+cat("\nAll checks passed.\n")
diff --git a/dev/benchmarks/frac-quart/verify_conc_disc.R b/dev/benchmarks/frac-quart/verify_conc_disc.R
new file mode 100644
index 000000000..c60fd14c9
--- /dev/null
+++ b/dev/benchmarks/frac-quart/verify_conc_disc.R
@@ -0,0 +1,73 @@
+# Brute-force check of concordant/discordant/independent-quartet formulae
+# for two binary bipartitions with 2x2 taxon table {p,q,r,s}.
+choose2 <- function(x) x * (x - 1) / 2
+
+# quartet resolution induced by a 0/1 labelling of 4 taxa: returns the
+# unordered pair-of-pairs if 2-2, else NA
+res <- function(lab) {
+  if (sum(lab) != 2) return(NA_character_)   # not 2-2
+  paste(sort(which(lab == lab[order(lab)][1]))[1:2], collapse = "")  # placeholder
+}
+# cleaner: canonical partition key for a 2-2 split of positions 1:4
+key22 <- function(side) {
+  if (length(unique(side)) != 2 || sum(side == side[1]) != 2) return(NA_character_)
+  a <- sort(which(side == side[1]))
+  paste(a, collapse = "")   # the pair on 'first' side identifies the 2|2 split up to complement
+}
+
+set.seed(1)
+ok <- TRUE
+for (rep in 1:2000) {
+  t <- sample(6:11, 1)
+  charState <- sample(0:1, t, replace = TRUE)   # I = state1
+  splitSide <- sample(0:1, t, replace = TRUE)   # A = side1
+  p <- sum(charState == 1 & splitSide == 1)
+  q <- sum(charState == 1 & splitSide == 0)
+  r <- sum(charState == 0 & splitSide == 1)
+  s <- sum(charState == 0 & splitSide == 0)
+
+  conc <- 0; disc <- 0; decisive <- 0
+  combs <- combn(t, 4)
+  for (j in seq_len(ncol(combs))) {
+    idx <- combs[, j]
+    ck <- key22(charState[idx]); sk <- key22(splitSide[idx])
+    if (!is.na(ck)) decisive <- decisive + 1
+    if (!is.na(ck) && !is.na(sk)) {
+      # same 2|2 split iff the pair-key matches OR is complementary
+      same <- (ck == sk) || (ck == paste(setdiff(1:4, as.integer(strsplit(ck,"")[[1]])), collapse=""))
+      if (same) conc <- conc + 1 else disc <- disc + 1
+    }
+  }
+  fConc <- choose2(p)*choose2(s) + choose2(q)*choose2(r)
+  fDisc <- p*q*r*s
+  n <- p + q
+  fDecisive <- choose2(n) * choose2(t - n)
+  if (conc != fConc || disc != fDisc || decisive != fDecisive) {
+    cat("MISMATCH pqrs=", p,q,r,s,
+        " conc", conc, fConc, " disc", disc, fDisc,
+        " dec", decisive, fDecisive, "\n"); ok <- FALSE; break
+  }
+  # fractional (independent) agreement with floors
+  pf <- function(x) pmax(x - 1, 0)
+  indConc <- pf(p)*pf(s) + pf(q)*pf(r)
+  indDec  <- pf(n)*pf(t - n)           # char's own independent quartets
+  if (indDec > 0) {
+    FQ <- indConc / indDec
+    if (FQ < -1e-9 || FQ > 1 + 1e-9) { cat("FQ out of range", FQ, p,q,r,s, "\n"); ok <- FALSE; break }
+  }
+  # self-agreement: char == split  => q=r=0 => FQ==1
+  if (q == 0 && r == 0 && indDec > 0 && abs(FQ - 1) > 1e-9) {
+    cat("self-agree != 1", FQ, p,q,r,s, "\n"); ok <- FALSE; break
+  }
+}
+cat(if (ok) "ALL CHECKS PASS\n" else "FAILED\n")
+
+# ind_decisive - ind_concordant identity (advisor: = pr + qs - 1 without floors)
+set.seed(2); id_ok <- TRUE
+for (i in 1:5000) {
+  p<-sample(0:6,1);q<-sample(0:6,1);r<-sample(0:6,1);s<-sample(0:6,1)
+  n<-p+q; tt<-p+q+r+s
+  lhs <- (n-1)*(tt-n-1) - ((p-1)*(s-1) + (q-1)*(r-1))
+  if (lhs != p*r + q*s - 1) { cat("identity fail\n"); id_ok<-FALSE; break }
+}
+cat(if (id_ok) "IDENTITY (no-floor) ind_dec - ind_conc == pr+qs-1: CONFIRMED\n" else "IDENTITY FAIL\n")
diff --git a/dev/benchmarks/frac-quart/verify_direction_dep.R b/dev/benchmarks/frac-quart/verify_direction_dep.R
new file mode 100644
index 000000000..cf3fcad4d
--- /dev/null
+++ b/dev/benchmarks/frac-quart/verify_direction_dep.R
@@ -0,0 +1,47 @@
+source("gfrank.R")   # reuse gf2rank (prints table; ignore)
+# Represent the SAME decisive quartet set two ways:
+#  - by the CHARACTER's pairing (4-cycle in K_{n,l}, I x O)      -> Delta_char
+#  - by the SPLIT's pairing     (4-cycle in K_{m,t-m}, A x B)     -> Delta_split
+# Concordant quartets pair the same either way; discordant differ.
+delta_both <- function(p, q, r, s) {
+  n <- p + q; l <- r + s; m <- p + r; b <- q + s   # split sizes
+  # taxon coords: state (I=1..n grouped p then q), side; build 4 cells
+  # cell membership: give each taxon (charState in {I,O}, side in {A,B})
+  taxa <- rbind(
+    cbind(I=1, A=1, id=seq_len(p)),                       # p: I&A
+    cbind(I=1, A=0, id=seq_len(q)),                       # q: I&B
+    cbind(I=0, A=1, id=seq_len(r)),                       # r: O&A
+    cbind(I=0, A=0, id=seq_len(s)))                       # s: O&B
+  N <- nrow(taxa)
+  # index taxa within char-graph: I-vertices and O-vertices
+  Iidx <- which(taxa[,"I"]==1); Oidx <- which(taxa[,"I"]==0)
+  Aidx <- which(taxa[,"A"]==1); Bidx <- which(taxa[,"A"]==0)
+  vc <- match(seq_len(N), Iidx); vo <- match(seq_len(N), Oidx)  # char graph coord
+  va <- match(seq_len(N), Aidx); vb <- match(seq_len(N), Bidx)  # split graph coord
+  edgeC <- function(iA, iB) (iA-1)*length(Oidx) + iB           # char edge id
+  edgeS <- function(a, b) (a-1)*length(Bidx) + b               # split edge id
+  rowsC <- list(); rowsS <- list()
+  quads <- combn(N, 4)
+  for (j in seq_len(ncol(quads))) {
+    t4 <- quads[, j]
+    isI <- taxa[t4,"I"]==1; isA <- taxa[t4,"A"]==1
+    if (sum(isI)!=2 || sum(isA)!=2) next                  # decisive = both resolve 2-2
+    Is <- t4[isI]; Os <- t4[!isI]; As <- t4[isA]; Bs <- t4[!isA]
+    # char pairing: (Is)|(Os)
+    ec <- c(edgeC(vc[Is[1]],vo[Os[1]]), edgeC(vc[Is[1]],vo[Os[2]]),
+            edgeC(vc[Is[2]],vo[Os[1]]), edgeC(vc[Is[2]],vo[Os[2]]))
+    es <- c(edgeS(va[As[1]],vb[Bs[1]]), edgeS(va[As[1]],vb[Bs[2]]),
+            edgeS(va[As[2]],vb[Bs[1]]), edgeS(va[As[2]],vb[Bs[2]]))
+    vC <- integer(length(Iidx)*length(Oidx)); vC[ec] <- 1; rowsC[[length(rowsC)+1]] <- vC
+    vS <- integer(length(Aidx)*length(Bidx)); vS[es] <- 1; rowsS[[length(rowsS)+1]] <- vS
+  }
+  mk <- function(L, w) if (length(L)) do.call(rbind, L) else matrix(0,0,w)
+  c(Dchar = gf2rank(mk(rowsC, length(Iidx)*length(Oidx))),
+    Dsplit= gf2rank(mk(rowsS, length(Aidx)*length(Bidx))),
+    char_wt = (n-1)*(l-1), split_wt = (m-1)*(b-1))
+}
+for (cfg in list(c(3,1,2,2), c(2,3,4,1), c(5,2,3,4), c(4,1,2,3))) {
+  d <- delta_both(cfg[1],cfg[2],cfg[3],cfg[4])
+  cat(sprintf("p,q,r,s=%-10s  Delta_char=%d (charWt=%d)   Delta_split=%d (splitWt=%d)\n",
+      paste(cfg,collapse=","), d["Dchar"], d["char_wt"], d["Dsplit"], d["split_wt"]))
+}
diff --git a/dev/benchmarks/frac-quart/verify_trit_ranks.R b/dev/benchmarks/frac-quart/verify_trit_ranks.R
new file mode 100644
index 000000000..f1830671c
--- /dev/null
+++ b/dev/benchmarks/frac-quart/verify_trit_ranks.R
@@ -0,0 +1,59 @@
+# Independent-quartet counts as GF(2) ranks in the character's bipartite
+# cycle space. A character quartet {i1,i2 | o1,o2} = a 4-cycle of K_{n,l}
+# (n = I-taxa, l = O-taxa); entailment = cycle addition (Nelson-Ladiges).
+gf2rank <- function(M) {                       # rank over GF(2)
+  if (!nrow(M) || !ncol(M)) return(0L)
+  M <- M %% 2; r <- 0L; ncol <- ncol(M); row <- 1L
+  for (col in seq_len(ncol)) {
+    piv <- which(M[row:nrow(M), col] == 1)
+    if (!length(piv)) next
+    piv <- piv[1] + row - 1L
+    tmp <- M[row, ]; M[row, ] <- M[piv, ]; M[piv, ] <- tmp
+    others <- which(M[, col] == 1); others <- others[others != row]
+    for (o in others) M[o, ] <- (M[o, ] + M[row, ]) %% 2
+    row <- row + 1L; r <- r + 1L
+    if (row > nrow(M)) break
+  }
+  r
+}
+
+# taxa: I-vertices 1..n (in A iff <=p), O-vertices 1..l (in A iff <=r)
+quartet_vec <- function(i1, i2, o1, o2, l) {
+  v <- integer(0)  # edge id = (i-1)*l + o
+  e <- c((i1-1)*l+o1, (i1-1)*l+o2, (i2-1)*l+o1, (i2-1)*l+o2)
+  e
+}
+counts <- function(p, q, r, s) {
+  n <- p + q; l <- r + s; nEdge <- n * l
+  inA_I <- function(i) i <= p
+  inA_O <- function(o) o <= r
+  rows_conc <- list(); rows_disc <- list(); rows_dec <- list()
+  for (i1 in 1:(n-1)) for (i2 in (i1+1):n) for (o1 in 1:(l-1)) for (o2 in (o1+1):l) {
+    nA <- inA_I(i1) + inA_I(i2) + inA_O(o1) + inA_O(o2)
+    if (nA != 2) next                              # split doesn't resolve -> not decisive
+    e <- c((i1-1)*l+o1, (i1-1)*l+o2, (i2-1)*l+o1, (i2-1)*l+o2)
+    v <- integer(nEdge); v[e] <- 1
+    # concordant iff the two I-taxa are on the same split side
+    conc <- (inA_I(i1) == inA_I(i2))
+    rows_dec[[length(rows_dec)+1]] <- v
+    if (conc) rows_conc[[length(rows_conc)+1]] <- v
+    else       rows_disc[[length(rows_disc)+1]] <- v
+  }
+  mk <- function(L) if (length(L)) do.call(rbind, L) else matrix(0, 0, nEdge)
+  list(conc = gf2rank(mk(rows_conc)),
+       disc = gf2rank(mk(rows_disc)),
+       dec  = gf2rank(mk(rows_dec)),
+       n_conc = length(rows_conc), n_disc = length(rows_disc))
+}
+pf <- function(x) pmax(x - 1, 0)
+cat(sprintf("%-14s %5s %5s %5s | %-22s %-10s\n",
+            "p,q,r,s", "conc", "disc", "dec", "A=(p1)(s1)+(q1)(r1)", "n-1 l-1"))
+for (cfg in list(c(2,1,1,2), c(3,2,2,3), c(3,1,2,2), c(2,2,2,2),
+                 c(4,1,1,4), c(3,3,1,1), c(2,3,4,1), c(5,2,3,4),
+                 c(1,2,3,4), c(4,0,0,4), c(3,0,2,3), c(2,1,0,3))) {
+  p<-cfg[1];q<-cfg[2];r<-cfg[3];s<-cfg[4]; n<-p+q; l<-r+s
+  x <- counts(p,q,r,s)
+  A <- pf(p)*pf(s) + pf(q)*pf(r)
+  cat(sprintf("%-14s %5d %5d %5d | A=%-3d disc? dec vs %-4d| char=(%d)(%d)=%d\n",
+      paste(cfg,collapse=","), x$conc, x$disc, x$dec, A, A + x$disc, n-1, l-1, (n-1)*(l-1)))
+}
diff --git a/dev/plans/frac-quart-weight-agreement.md b/dev/plans/frac-quart-weight-agreement.md
new file mode 100644
index 000000000..7fe19b3d4
--- /dev/null
+++ b/dev/plans/frac-quart-weight-agreement.md
@@ -0,0 +1,272 @@
+# Fractional (Nelson–Ladiges) quartet weighting for `QuartetConcordance()`
+
+Branch: `frac-quart`. Goal: port Nelson & Ladiges (1992, *Syst. Biol.* 41:490–494)
+fractional weighting of three-item statements into the **unrooted quartet**
+concordance measure, so that logically redundant quartets are counted as the
+information they actually carry (trits) rather than as raw combinatorial volume.
+
+This file is the settled design. Everything below was derived and, where noted,
+verified by brute-force enumeration / GF(2) rank (scratch scripts, this session).
+
+---
+
+## 1. Decisions locked in
+
+- **Unrooted quartets only.** No polarity / no "informative state" asymmetry.
+  N&L's rooted `2/n` symmetrises to the quartet weight below.
+- **New arg `unit = c("quartet", "trit")`** on `QuartetConcordance()`. It
+  *replaces* any `fractional` flag. `"quartet"` = current behaviour (raw counts);
+  `"trit"` = N&L fractional (independent-quartet) currency.
+- **Normalization = coverage by the *reported* unit** (option "b"): only a
+  character-split **identical** to the tree split scores full marks; nesting and
+  crossing get partial credit. See §4.
+- **`weight` (TRUE/FALSE) is retained and orthogonal to `unit`.** It controls the
+  *pooling amount*, not the per-pair quality. See §5.
+- **`return` (edge/char/tree) retained;** it selects which unit is "reported" and
+  therefore which weight normalises. See §4–5.
+- This slots into the package's existing **quality × amount** decomposition
+  (`ConcordanceTable()` / `QACol()`): quality = per-pair coverage ratio, amount =
+  N&L fractional weight (trit content).
+
+---
+
+## 2. Currency: trits, quartets = 4-cycles
+
+For `t` taxa, a binary character (no missing data) is a bipartition of sizes
+`(k, t−k)`. The quartets it resolves are the **4-cycles of the complete bipartite
+graph `K_{k, t−k}`** (two taxa each side = one 4-cycle). Nelson–Ladiges
+entailment `(ab,cd) + (ab,ce) → (ab,de)` is exactly GF(2) cycle addition, so the
+independent-quartet count is the cyclomatic number:
+
+```
+independent quartets of a split of sizes (k, t−k)  =  (k−1)(t−k−1)      # "trits"
+raw quartets it resolves                           =  C(k,2)·C(t−k,2)
+per-quartet fractional weight                      =  4 / (k(t−k))       # = indep/raw
+```
+
+`(k−1)(t−k−1)` is the "fractional weight" `W` of the whole character/split. It is
+`0` for a constant (`k=0`) or autapomorphy (`k=1`) — so **uninformative characters
+self-zero in trit currency**, no special-casing needed.
+
+Smallest N&L case (t=5, k=2): `4/(2·3)·3 = 2` trits from 3 raw quartets. ✓
+
+**Multistate (solved, not deferred).** A quartet is decisive only when two taxa
+share one state and two share another, so every decisive quartet lives in the
+`K_{n_i, n_j}` edge-block between two states `i, j`. Those blocks are
+edge-disjoint, and the N&L entailment never crosses them (it forces the shared
+`c,d,e` taxa into one state), so **trits add over state pairs**:
+
+```
+W_char = Σ_{i<j} (n_i−1)(n_j−1)      A = Σ_{i<j} A_ij      weight_{ij} = 4/(n_i n_j)
+```
+
+reducing to `(n−1)(t−n−1)` for binary (one pair). Each state-pair is an
+independent binary sub-problem, so the binary machinery below is applied
+per state-pair and summed; a binary character is the single-pair case. This is
+the *same* trit currency (multistate and binary are directly comparable), not a
+separate unit. Verified by GF(2) rank of the decisive **and** concordant quartet
+sets on random 2/3/4-state characters (`dev/benchmarks/frac-quart/multistate_trit.R`).
+
+---
+
+## 3. Per-pair structure: the 2×2 table
+
+Character `c` (bipartition, `n = n_c` taxa in one state) vs tree split `k`
+(sizes `m = m_k`, `t−m`). Cross-tabulate into cells:
+
+```
+            split A    split B
+state I       p          q         (n   = p+q)
+state O       r          s         (t−n = r+s)
+                                    (m   = p+r,  t−m = q+s)
+```
+
+**Kernel decode** (verified against `src/quartet_concordance.cpp`, binary map
+state1=I, side1=A ⇒ `n1[1]=p, n0[1]=q, n1[0]=r, n0[0]=s`):
+
+```
+raw concordant  conc(k,c) = C(p,2)C(s,2) + C(q,2)C(r,2)
+raw decisive    dec(k,c)  = conc + p·q·r·s              # = concordant + discordant (both-resolve)
+```
+
+**Trit (independent) quantities** — GF(2) rank of the 4-cycle sets, verified by
+enumeration:
+
+```
+concordant trits   A(k,c) = (p−1)₊(s−1)₊ + (q−1)₊(r−1)₊      # ₊ = max(·,0); floors are load-bearing
+```
+
+The floors matter: without them the identity `ind_dec − ind_conc = pr + qs − 1`
+pushes self-agreement above 1.
+
+**Direction-dependence of the decisive count** (verified): the decisive set's
+trit rank differs by which partition supplies the pairing, because discordant
+quartets are paired differently by character vs split —
+`Δ_char = (n−1)(t−n−1)`, `Δ_split = (m−1)(t−m−1)` (generic, all four cells ≥ 1).
+`A` is symmetric (shared). This is *why* normalization depends on `return`.
+
+**Three regimes by number of empty cells** (standard split-compatibility trichotomy):
+
+| empty cells | relationship | discordant | verdict |
+|---|---|---|---|
+| 2 | **identical** (char split = tree split) | none | full support |
+| 1 | **nested** (compatible, proper subset) | none | partial support (real!) |
+| 0 | **crossing** (incompatible) | `pqrs > 0` | genuine test; conflict lowers score |
+
+Nesting is a *passed test*, not homoplasy: no evidence against ⇒ evidence for
+(cf. viviparity → Mammalia despite the platypus). It must get partial positive
+credit — not exclusion.
+
+---
+
+## 4. The measure (option b): coverage by the reported unit
+
+Per pair, quality = concordant trits over the **reported unit's** own weight:
+
+```
+edge  (report split k):    Q(k,c) = A(k,c) / W_k,   W_k = (m−1)(t−m−1)
+char  (report char  c):    Q(k,c) = A(k,c) / W_c,   W_c = (n−1)(t−n−1)
+```
+
+`A ≤ W_reported` always (⇒ `Q ∈ [0,1]`), and **per state-pair**
+`Q_ij = 1 ⟺ A_ij = W_reported ⟺ identical sub-bipartition`. So for a **binary**
+character, only an identical character-split scores full marks; nesting and
+crossing are partial, with conflict strictly lowering the score. Verified ladder
+(t=8, edge):
+
+| (p,q,r,s) | relationship | A | W_k | Q |
+|---|---|---|---|---|
+| 4,0,0,4 | identical | 9 | 9 | 1.00 |
+| 3,0,2,3 | nested | 4 | 8 | 0.50 |
+| 3,1,1,3 | mild crossing | 4 | 9 | 0.44 |
+| 2,2,2,2 | max crossing | 2 | 9 | 0.22 |
+
+**Multistate: `Q = 1 ⟺ the split is displayed by the character`** (§5 pools the
+per-pair `Q_ij` by `M = min(W_c, W_k)`). When every state block lies wholly on
+one side of the split, each split-relevant pair scores `Q_ij = 1` and the one
+split-orthogonal pair has `W_k^{ij} = 0 ⇒ M = 0` and drops out, so the pooled
+`Q = 1` — even though `A_tot < W_c,tot` and the multistate character is *not*
+"identical" to the (binary) split. This is the correct generalisation of
+"identical" (for a binary character, displayed = identical), and the intended
+behaviour: a reproductive character coded `{oviparous | viviparous | ovovivip.}`
+that cleanly refines Mammalia fully supports the Mammalia split. A state block
+that *straddles* the split makes its pairs crossing (`Q_ij < 1`), lowering the
+pooled score. Verified: `t7,t8,t9`-refining char scores `1.0` on every displayed
+split, `0.667` on splits that cut a block (see `test-Concordance.R`).
+
+Raw-currency (`unit="quartet"`) analogue: `Q = conc / (raw quartets of the
+reported unit)` — same shape, `conc` and `C(m,2)C(t−m,2)` instead of `A`, `W_k`.
+
+---
+
+## 5. Pooling (`weight` acts on the amount, not the quality)
+
+`W_k` is constant across characters for a fixed edge, so `weight` must act on the
+**amount** (pooling weight), else it goes inert:
+
+```
+edge FQ(k) = Σ_c  M(k,c)·Q(k,c)  /  Σ_c  M(k,c)
+char FQ(c) = Σ_k  M(k,c)·Q(k,c)  /  Σ_k  M(k,c)
+```
+
+- `M(k,c)` = shared information = **`min(W_c, W_k)`** (the hBest analogue used by
+  `ClusteringConcordance`). Preserves "only identical = full marks" (ceiling needs
+  `A = W_c = W_k`).
+- `weight = FALSE` ⇒ `M ≡ 1` (per-pair-equal mean of `Q`).
+- `unit = "quartet"` ⇒ `M` = raw shared-quartet count; `Q` uses raw conc.
+
+**Quality × amount map to the existing table:** `Q` is `QACol`'s *quality* axis;
+`W_c` (trit content) is the *amount* axis. The trit port reuses
+`ConcordanceTable()` rather than parallelling it.
+
+---
+
+## 6. Build plan
+
+1. **R prototype first** (binary, no missing). Compute cells `{p,q,r,s}` per
+   (split, char) from the split logical vectors + character integer vectors;
+   derive `A`, `W_c`, `W_k`, `Q`, `M`; pool per §5. Keep the existing C++ kernel
+   for the raw path; the trit path can start in R since cells are a cheap crosstab.
+2. **Validate**: (i) brute-force check `A` and the ladder on toy cells;
+   (ii) run on `congreveLamsdellMatrices` — report edge-score movement
+   `quartet → trit`, the identical/nested/crossing census per edge, and a
+   quality/amount table; (iii) confirm `unit="quartet"` reproduces current
+   `QuartetConcordance` bit-for-bit.
+3. **Kernel** (later, perf only): the R per-state-pair path is already correct
+   and validated for binary + multistate + missing data, so a C++ port of `A`
+   is now a *performance* optimisation (large matrices), not a correctness
+   prerequisite.
+4. **Tests**: extend `tests/testthat/test-QuartetConcordance*` — `unit` switch,
+   identical→1, nested partial, crossing < nested, uninformative→dropped,
+   `unit="quartet"` == legacy.
+
+## 7. Deferred / open
+
+- **Multistate**: ~~deferred~~ **DONE.** Trits add over state pairs
+  (`W_char = Σ_{i<j}(n_i−1)(n_j−1)`, `A = Σ A_ij`, per-pair weight `4/(n_i n_j)`);
+  see §2. Implemented as a per-state-pair loop in `.TritConcordance()`
+  (`R/Concordance.R`); GF(2)-verified; binary is the single-pair special case.
+- **Missing data**: ~~deferred~~ handled by construction — the per-state-pair
+  cell counts restrict to scored taxa (`scored <- !is.na(col)`), giving each
+  character its own effective `t`, `W_c` and per-pair `W_k`. Validated on a
+  toy with `?`/`-`/ambiguous tokens (`validate_trit.R`), but the validation
+  corpus (congreveLamsdell) has no missing data, so the semantics are proven
+  self-consistent (`Q∈[0,1]`, pkg==independent ref) rather than benchmarked.
+- **Random-expectation baseline**: ~~deferred~~ **DONE** (`unit = "trit"`).
+  New arg `normalize = FALSE / TRUE / <int>` on `QuartetConcordance()`, mirroring
+  `ClusteringConcordance()`. The crossing floor (≈ 0.22 in §4) is not zero, so
+  without correction a maximally conflicting character still scores positive and
+  the floor is split-size dependent; `normalize` re-zeros against the
+  concordance *expected by chance*, exactly as `ClusteringConcordance()` re-zeros
+  MI against `miRand`.
+  - **Null model (chosen): fixed-marginal randomization** — reassign each
+    character's tokens across the leaves at random, holding its state counts and
+    the split sizes fixed (the multivariate-hypergeometric confusion table; the
+    same null as `ClusteringConcordance`'s `miRand` and `Consistency`'s
+    `ExpectedLength` tip-shuffle). This is the *only* coherent null for the trit
+    currency: it re-zeros each `(split, char, state-pair)` against its own
+    combinatorial floor and extends unchanged to multistate and missing data.
+  - **Rejected alternative: the flat Minh-style `1/3`** (three quartet
+    resolutions ⇒ expected concordance `1/3`). It is a *raw-quartet* heuristic
+    that ignores marginals and does **not** extend to floored trit counts (the
+    trit floor is `0.22`, not `1/3`), so it cannot re-zero the trit measure
+    coherently. Its only remaining use would be re-baselining the *raw* `quartet`
+    unit — see the raw-unit note below.
+  - **Two estimators** (mirroring `ClusteringConcordance`): `normalize = TRUE`
+    computes the **exact** expected pools from the hypergeometric pmf (a small
+    double sum, cached on `(n_i, n_j, M, t)`); `normalize = <int>` estimates them
+    by **Monte-Carlo** tip-shuffle. Validated to agree within MC error
+    (`dev/benchmarks/frac-quart/chance_baseline.R` at the cell level;
+    `test-Concordance.R` end-to-end).
+  - **What varies under the null.** Only `A` (and, for *multistate* pairs where a
+    pair's side-A count `mA` is itself random, `wk` and `M = min(w_c, w_k)`) vary;
+    for *binary* characters `mA = M`, `tP = t` so `wk` is fixed and only `A`
+    varies (the clean `ClusteringConcordance` case). The exact estimator therefore
+    accumulates `E[m]`, `E[m·A/w_k]`, `E[m·A/w_c]` and re-zeros the *pooled*
+    edge/char ratio against the pooled expectation, with the **same** `weight`
+    aggregation as the observed score.
+  - **Guarantees / guards.** `.Rezero(1, z) = 1`, so a **displayed** split still
+    scores exactly `1` (the ceiling invariant of §4 survives chance correction —
+    only the floor is lifted); conflicting characters go **negative** ("below
+    random", as in `ClusteringConcordance`). Values are returned **unclamped**
+    (may fall below `−1`); clamp to `[−1, 1]` only at plot time (`QCol`/`QACol`),
+    never inside the measure. The `z → 1` blow-up is guarded (→ `NA`).
+    `normalize` defaults to `FALSE`, so the published measure is unchanged unless
+    the user opts in.
+  - **Raw `quartet` unit: DONE.** `normalize` now works for `unit = "quartet"`
+    too, using the *same* fixed-marginal null. Per state-pair `conc` and `dec`
+    are polynomials in the cells (no floors), so the exact expectation
+    `E[conc]`, `E[dec]` is a clean sum over the same trivariate-hypergeometric
+    pmf as `.ExpectedTrit` (`.ExpectedQuartet` / `.QuartetExpect`); MC via
+    `.QuartetMC` reshuffles tokens and re-scores through the kernel. The pooled
+    `conc/dec` ratio is re-zeroed against `E[conc]/E[dec]` (weighted) or with
+    the same cell-matching as trits (unweighted). Validated exact≈MC across
+    edge/char × weight × binary/multistate (`test-Concordance.R`), and
+    `normalize = FALSE` is byte-identical to the published raw measure. The flat
+    `1/3` baseline stays rejected (it does not adapt to marginals; the
+    fixed-marginal null generalises it).
+- **`M = min(W_c,W_k)`** (settled): the pooling amount is the *shared*
+  information, so it must be symmetric between character and split (both
+  `return`s pool by the same `M`); `M = W_k` would make the char return weight a
+  pair by the *edge's* content, which is incoherent as shared information.
+  `min(W_c,W_k)` also drops into `ConcordanceTable`'s `hBest·n` slot (§5).
diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib
index 22e2babfb..8fd1e150c 100644
--- a/inst/REFERENCES.bib
+++ b/inst/REFERENCES.bib
@@ -1,3 +1,5 @@
+@comment{List entries alphabetically by key}
+
 @article{Aberer2013,
   title = {Pruning rogue taxa improves phylogenetic accuracy: an efficient algorithm and webservice},
   author = {Aberer, Andre J. and Krompass, Desmithconsnis and Stamatakis, Alexandros},
@@ -506,6 +508,17 @@ @article{Muller2005
   number = {1}
 }
 
+@article{Nelson1992,
+  title = {Information content and fractional weight of three-taxon statements},
+  author = {Nelson, Gareth and Ladiges, Pauline Y.},
+  year = {1992},
+  journal = {Systematic Biology},
+  volume = {41},
+  number = {4},
+  pages = {490--494},
+  doi = {10.1093/sysbio/41.4.490}
+}
+
 @article{Nguyen2015,
   title = {{IQ-TREE}: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies},
   author = {Nguyen, Lam-Tung and Schmidt, Heiko A. and von Haeseler, Arndt and Minh, Bui Quang},
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 8d0dc270b..733e626c7 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -179,6 +179,7 @@ TreeDist
 TreeTools
 Trueman
 UHL
+unclamped
 VERBOOM
 VILHELMSEN
 VILLAR
diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd
index 9d6a0831b..61b010f79 100644
--- a/man/ConcordanceTable.Rd
+++ b/man/ConcordanceTable.Rd
@@ -11,7 +11,7 @@ ConcordanceTable(
   largeClade = 0,
   xlab = "Edge",
   ylab = "Character",
-  normalize = TRUE,
+  chanceCorrect = TRUE,
   plot = TRUE,
   marginSize = 0L,
   paintSize = 0L,
@@ -36,17 +36,20 @@ at edges whose descendants are both contain more than \code{largeClade} leaves.}
 
 \item{ylab}{Character giving a label for the y axis.}
 
-\item{normalize}{Controls how the \emph{expected} mutual information (the zero
-point of the scale) is determined.
+\item{chanceCorrect}{Sets the zero point of the scale.
 \itemize{
-\item \code{FALSE}: no chance correction; MI is scaled only by its maximum.
-\item \code{TRUE}: subtract the analytical expected MI for random association.
-\item \verb{<integer>}: subtract an empirical expected MI estimated from that
-number of random trees.
-}
-
-In all cases, 1 corresponds to the maximal attainable MI for the pair
-(\code{hBest}), and 0 corresponds to the chosen expectation.}
+\item If \code{FALSE}, zero corresponds to zero MI.
+\item If \code{TRUE}, zero is the value expected when each character's tokens are
+reassigned at random across the leaves, holding its state frequencies and the
+split sizes fixed: \code{QuartetConcordance()} computes this expectation exactly,
+whereas \code{ClusteringConcordance()} approximates it, accurately for large trees
+(~200+ taxa) but neglecting correlation between splits.
+\item If a positive integer \code{n}, the expectation is sampled over \code{n} random
+reassignments of each character's tokens (\code{QuartetConcordance()}), or against
+\code{n} uniformly random trees (\code{ClusteringConcordance()}).
+\emph{Hint: Clamp chance-corrected values to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} /
+\code{\link[=QACol]{QACol()}}}.
+}}
 
 \item{plot}{Logical specifying whether to draw the plot.}
 
diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd
index 4436e882d..c29a9880c 100644
--- a/man/SiteConcordance.Rd
+++ b/man/SiteConcordance.Rd
@@ -9,11 +9,18 @@
 \alias{SharedPhylogeneticConcordance}
 \title{Concordance factors}
 \usage{
-ClusteringConcordance(tree, dataset, return = "edge", normalize = TRUE)
+ClusteringConcordance(tree, dataset, return = "edge", chanceCorrect = TRUE)
 
 MutualClusteringConcordance(tree, dataset)
 
-QuartetConcordance(tree, dataset = NULL, weight = TRUE, return = "edge")
+QuartetConcordance(
+  tree,
+  dataset = NULL,
+  weight = TRUE,
+  return = "edge",
+  unit = c("nrqs", "quartet"),
+  chanceCorrect = TRUE
+)
 
 PhylogeneticConcordance(tree, dataset)
 
@@ -33,24 +40,35 @@ by each character's information content;
 \item \code{"tree"}: an overall tree‑level concordance score;
 \item \code{"all"}: a full array of MI components and normalized values for every
 split–character pair.
-}
-
-Matching is case‑insensitive and partial.}
+}}
 
-\item{normalize}{Controls how the \emph{expected} mutual information (the zero
-point of the scale) is determined.
+\item{chanceCorrect}{Sets the zero point of the scale.
 \itemize{
-\item \code{FALSE}: no chance correction; MI is scaled only by its maximum.
-\item \code{TRUE}: subtract the analytical expected MI for random association.
-\item \verb{<integer>}: subtract an empirical expected MI estimated from that
-number of random trees.
-}
-
-In all cases, 1 corresponds to the maximal attainable MI for the pair
-(\code{hBest}), and 0 corresponds to the chosen expectation.}
+\item If \code{FALSE}, zero corresponds to zero MI.
+\item If \code{TRUE}, zero is the value expected when each character's tokens are
+reassigned at random across the leaves, holding its state frequencies and the
+split sizes fixed: \code{QuartetConcordance()} computes this expectation exactly,
+whereas \code{ClusteringConcordance()} approximates it, accurately for large trees
+(~200+ taxa) but neglecting correlation between splits.
+\item If a positive integer \code{n}, the expectation is sampled over \code{n} random
+reassignments of each character's tokens (\code{QuartetConcordance()}), or against
+\code{n} uniformly random trees (\code{ClusteringConcordance()}).
+\emph{Hint: Clamp chance-corrected values to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} /
+\code{\link[=QACol]{QACol()}}}.
+}}
 
 \item{weight}{Logical specifying whether to weight sites according to the
-number of quartets they are decisive for.}
+quartet content that they share with each split.}
+
+\item{unit}{Character specifying the currency in which quartets are counted:
+\itemize{
+\item \code{"nrqs"} (default): non-redundant quartet statements.
+At a given split, only a character whose own split is identical scores
+full marks; nested (compatible) characters receive partial support, and
+incompatible characters score lower still.
+\item \code{"quartet"}: each resolved quartet counts once,
+analogous to the site concordance factor \insertCite{Minh2020}{TreeSearch}.
+}}
 }
 \value{
 \code{ClusteringConcordance(return = "all")} returns a 3D array where each
@@ -71,7 +89,7 @@ expectation. \code{NA} is returned when \code{hBest = 0} (no information potenti
 corresponds to a split (an edge of the tree) and gives the normalized mutual
 information between that split and the character data, averaged across all
 characters.
-When \code{normalize = TRUE} (default), values are scaled relative to random
+When \code{chanceCorrect = TRUE} (default), values are scaled relative to random
 expectation; when \code{FALSE}, raw mutual information normalized by \code{hBest} is
 returned.
 
@@ -103,6 +121,10 @@ concordance index calculated at each site, averaged across all splits.
 Unlike the \code{"edge"} result, this vector is unnamed: characters retain no
 persistent identifier once collapsed to patterns.
 
+With \code{unit = "nrqs"}, per-character values are much smaller than 1 even
+where a character perfectly matches a split in a tree. Values rank characters
+on a given tree, but are not comparable between trees.
+
 \code{PhylogeneticConcordance()} returns a numeric vector giving the
 phylogenetic information of each split in \code{tree}, named according to the
 split's internal numbering.
@@ -130,22 +152,6 @@ is quantified using mutual information (MI). All reported values are scaled
 so that 1 corresponds to the maximum possible mutual information for each
 split–character pair (\code{hBest}).
 
-The \code{normalize} argument specifies how the zero point is defined.
-\itemize{
-\item If \code{normalize = FALSE}, zero corresponds to \emph{zero} MI, without correcting
-for the positive bias that arises because MI is rarely exactly zero in
-finite samples.
-\item If \code{normalize = TRUE}, the expected MI is computed using an analytical
-approximation based on the distribution of character tokens. This is fast
-and generally accurate for large trees (~200+ taxa), but does not account
-for correlation between splits.
-\item If \code{normalize} is a positive integer \code{n}, the expected MI is estimated
-empirically by fitting each character to \code{n} uniformly random trees and
-averaging the resulting MI values. This Monte Carlo approach provides a
-more accurate baseline for small trees, for which the analytical
-approximation is biased. Monte Carlo standard errors are returned.
-}
-
 \code{MutualClusteringConcordance()} provides a character‑wise summary that
 emphasises each character’s best‑matching split(s). It treats each character
 as a simple tree and computes the mutual clustering information between this
@@ -153,9 +159,9 @@ character‑tree and the supplied phylogeny. High values identify characters
 whose signal is well represented anywhere in the tree, even if concentrated
 on a single edge.
 
-\code{QuartetConcordance()} is the proportion of quartets (sets of four leaves)
-that are decisive for a split which are also concordant with it
-For example, a quartet with the characters \verb{0 0 0 1} is not decisive, as
+\code{QuartetConcordance()} measures the agreement between each character and each
+split in the currency of quartets (sets of four leaves).
+A quartet with the characters \verb{0 0 0 1} is not decisive, as
 all relationships between those leaves are equally parsimonious.
 But a quartet with characters \verb{0 0 1 1} is decisive, and is concordant
 with any tree that groups the first two leaves together to the exclusion
@@ -165,12 +171,21 @@ In contrast to the site concordance factor
 quartets that are decisive for a branch.
 Doing so circumvents the criticisms of \insertCite{Goloboff2024;textual}{TreeSearch}.
 
-By default, the reported value weights each site by the number of quartets
-it is decisive for.  This value can be interpreted as the proportion of
-all decisive quartets that are concordant with a split.
+The quartets that a split resolves are not logically independent: because
+\eqn{(ab, cd)} and \eqn{(ab, ce)} jointly entail \eqn{(ab, de)}
+\insertCite{Nelson1992}{TreeSearch}, a split of sizes \eqn{(k, t - k)}
+resolves \eqn{\binom{k}{2}\binom{t - k}{2}} quartets, of which only
+\eqn{(k - 1)(t - k - 1)} are non-redundant.
+\code{unit} selects which of the two is counted.  By default, agreement is
+measured against the non-redundant content, so that a character scores full
+marks only where it displays the split, rather than merely agreeing with a
+large combinatorial volume of its quartets.
+
+By default, the reported value weights each site by the quartet content it
+shares with the split.
 If \code{weight = FALSE}, the reported value is the mean of the concordance
 value for each site.
-Consider a split associated with two sites:
+Consider a split associated with two sites (counting in \code{unit = "quartet"}):
 one that is concordant with 25\% of 96 decisive quartets, and
 a second that is concordant with 75\% of 4 decisive quartets.
 If \code{weight = TRUE}, the split concordance will be 24 + 3 / 96 + 4 = 27\%.
@@ -247,6 +262,8 @@ myMatrix <- congreveLamsdellMatrices[[10]]
 ClusteringConcordance(TreeTools::NJTree(myMatrix), myMatrix)
 }
 \references{
+\insertAllCited{}
+
 \insertAllCited{}
 }
 \seealso{
diff --git a/man/expected_mi.Rd b/man/expected_mi.Rd
index 9b455127e..21159a2cd 100644
--- a/man/expected_mi.Rd
+++ b/man/expected_mi.Rd
@@ -22,7 +22,7 @@ partitions of the same \code{N} items, under the hypergeometric null in which th
 block sizes (marginals) of each partition are fixed but the items are
 associated at random.  Subtracting this baseline from an observed mutual
 information yields a chance-corrected ("adjusted") mutual information, as
-applied by \code{\link{SiteConcordance}}\code{(normalize = TRUE)}; it is most material for
+applied by \code{\link{SiteConcordance}}\code{(chanceCorrect = TRUE)}; it is most material for
 small or unbalanced partitions, where raw mutual information is appreciably
 inflated by chance agreement.
 }
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index fb9bb0149..3d1f0ca32 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -31,6 +31,30 @@ BEGIN_RCPP
     return R_NilValue;
 END_RCPP
 }
+// quartet_expect
+List quartet_expect(const LogicalMatrix splits, const IntegerMatrix characters);
+RcppExport SEXP _TreeSearch_quartet_expect(SEXP splitsSEXP, SEXP charactersSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const LogicalMatrix >::type splits(splitsSEXP);
+    Rcpp::traits::input_parameter< const IntegerMatrix >::type characters(charactersSEXP);
+    rcpp_result_gen = Rcpp::wrap(quartet_expect(splits, characters));
+    return rcpp_result_gen;
+END_RCPP
+}
+// nrqs_expect
+List nrqs_expect(const LogicalMatrix splits, const IntegerMatrix characters);
+RcppExport SEXP _TreeSearch_nrqs_expect(SEXP splitsSEXP, SEXP charactersSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const LogicalMatrix >::type splits(splitsSEXP);
+    Rcpp::traits::input_parameter< const IntegerMatrix >::type characters(charactersSEXP);
+    rcpp_result_gen = Rcpp::wrap(nrqs_expect(splits, characters));
+    return rcpp_result_gen;
+END_RCPP
+}
 // expected_mi
 double expected_mi(const IntegerVector& ni, const IntegerVector& nj);
 RcppExport SEXP _TreeSearch_expected_mi(SEXP niSEXP, SEXP njSEXP) {
diff --git a/src/TreeSearch-init.c b/src/TreeSearch-init.c
index 66903e2d7..c7ecc6e12 100644
--- a/src/TreeSearch-init.c
+++ b/src/TreeSearch-init.c
@@ -21,6 +21,8 @@ extern SEXP _TreeSearch_mi_key(SEXP, SEXP);
 
 // extern SEXP _TreeSearch_astar_search_r(SEXP, SEXP, SEXP);
 extern SEXP _TreeSearch_quartet_concordance(SEXP, SEXP);
+extern SEXP _TreeSearch_quartet_expect(SEXP, SEXP);
+extern SEXP _TreeSearch_nrqs_expect(SEXP, SEXP);
 extern SEXP _TreeSearch_ts_fitch_score(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP _TreeSearch_ts_na_char_steps(SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP _TreeSearch_ts_char_steps(SEXP, SEXP, SEXP, SEXP, SEXP);
@@ -76,6 +78,8 @@ static const R_CallMethodDef callMethods[] = {
 
   // {"_TreeSearch_astar_search_r",    (DL_FUNC) &_TreeSearch_astar_search_r, 3},
   {"_TreeSearch_quartet_concordance",(DL_FUNC) &_TreeSearch_quartet_concordance, 2},
+  {"_TreeSearch_quartet_expect",    (DL_FUNC) &_TreeSearch_quartet_expect, 2},
+  {"_TreeSearch_nrqs_expect",       (DL_FUNC) &_TreeSearch_nrqs_expect, 2},
   {"_TreeSearch_ts_fitch_score",    (DL_FUNC) &_TreeSearch_ts_fitch_score, 12},
   {"_TreeSearch_ts_na_char_steps", (DL_FUNC) &_TreeSearch_ts_na_char_steps, 5},
   {"_TreeSearch_ts_char_steps", (DL_FUNC) &_TreeSearch_ts_char_steps, 5},
diff --git a/src/concordance_expect.cpp b/src/concordance_expect.cpp
new file mode 100644
index 000000000..b6b201920
--- /dev/null
+++ b/src/concordance_expect.cpp
@@ -0,0 +1,215 @@
+#include <Rcpp.h>
+#include <vector>
+#include <algorithm>
+#include <cmath>
+
+using namespace Rcpp;
+
+// Expected concordance under the fixed-marginal (hypergeometric) null used by
+// QuartetConcordance(chanceCorrect = TRUE).  For one character each token is
+// reassigned at random across the scored leaves, holding the state counts and
+// the split sizes fixed.  For a state-pair (i, j) with counts nI, nJ over t
+// scored leaves, of which M are on side A of a split, the 2x2 cell vector
+// (p, q, r, s) is trivariate-hypergeometric:
+//   P(p, r) = C(nI,p) C(nJ,r) C(t-nI-nJ, M-p-r) / C(t, M)
+// with q = nI - p, s = nJ - r, and mA = p + r on side A.  These kernels sum the
+// per-pair contributions exactly (double sum over p, r), mirroring the R
+// reference `.ExpectedQuartet`/`.QuartetExpect` and `.ExpectedTrit`/
+// `.CharTritExpect` in R/Concordance.R.  The p-outer / r-inner accumulation
+// order matches the R code so the results agree to machine precision.
+
+// C(z, 2) with the same z < 2 -> 0 guard as the R `choose2`.
+static inline double choose2(double z) {
+  return z < 2.0 ? 0.0 : z * (z - 1.0) / 2.0;
+}
+
+// (z - 1) floored at 0, matching the R `pos(z - 1)`.
+static inline double posm1(double z) {
+  return z < 1.0 ? 0.0 : z - 1.0;
+}
+
+// Per-character scored state counts + per-split side-A size, shared by both
+// kernels.  `stateOf[t]` is the compact 0-based state index of a scored leaf, or
+// -1 if the leaf is NA for this character.
+struct CharCounts {
+  int tc;                       // number of scored leaves
+  std::vector<int> cnt;         // count of each compact state
+  std::vector<int> mSideA;      // per split: scored leaves on side A
+};
+
+static CharCounts char_counts(const LogicalMatrix& splits,
+                              const IntegerMatrix& characters, int c) {
+  const int n_taxa = splits.nrow();
+  const int n_splits = splits.ncol();
+  CharCounts out;
+  out.tc = 0;
+  out.mSideA.assign(n_splits, 0);
+
+  // Map raw state values to compact 0-based indices (sorted ascending, as R's
+  // sort(unique(...)) does; the pair loop is symmetric so order only affects
+  // which is "i" vs "j", not the sums).
+  std::vector<int> raw;
+  raw.reserve(n_taxa);
+  std::vector<int> stateOf(n_taxa, -1);
+  for (int t = 0; t < n_taxa; ++t) {
+    int st = characters(t, c);
+    if (IntegerVector::is_na(st)) continue;
+    ++out.tc;
+    raw.push_back(st);
+  }
+  std::sort(raw.begin(), raw.end());
+  raw.erase(std::unique(raw.begin(), raw.end()), raw.end());
+  out.cnt.assign(raw.size(), 0);
+
+  for (int t = 0; t < n_taxa; ++t) {
+    int st = characters(t, c);
+    if (IntegerVector::is_na(st)) continue;
+    int idx = int(std::lower_bound(raw.begin(), raw.end(), st) - raw.begin());
+    stateOf[t] = idx;
+    ++out.cnt[idx];
+    for (int s = 0; s < n_splits; ++s) {
+      if (splits(t, s)) ++out.mSideA[s];
+    }
+  }
+  return out;
+}
+
+// Exact E[conc], E[dec] for one state-pair over all needed side-A sizes.
+// Returns, for each M in 0..tc, the pair's (Econc, Edec); callers pick M[split].
+static void expected_quartet_by_M(int nI, int nJ, int tc,
+                                  std::vector<double>& eConc,
+                                  std::vector<double>& eDec) {
+  const int nOther = tc - nI - nJ;
+  eConc.assign(tc + 1, 0.0);
+  eDec.assign(tc + 1, 0.0);
+  // Hoist the log-choose terms out of the M / p / r loops (the dominant cost).
+  std::vector<double> lcI(nI + 1), lcJ(nJ + 1), lcO(nOther + 1), lcT(tc + 1);
+  for (int p = 0; p <= nI; ++p) lcI[p] = R::lchoose(nI, p);
+  for (int r = 0; r <= nJ; ++r) lcJ[r] = R::lchoose(nJ, r);
+  for (int o = 0; o <= nOther; ++o) lcO[o] = R::lchoose(nOther, o);
+  for (int M = 0; M <= tc; ++M) lcT[M] = R::lchoose(tc, M);
+  for (int M = 0; M <= tc; ++M) {
+    const double lt = lcT[M];
+    double aConc = 0.0, aDec = 0.0;
+    const int pMax = std::min(nI, M);
+    for (int p = 0; p <= pMax; ++p) {
+      const int q = nI - p;
+      const double c_p = choose2(p), c_q = choose2(q), lp = lcI[p];
+      const int rMax = std::min(nJ, M - p);
+      for (int r = 0; r <= rMax; ++r) {
+        const int oA = M - p - r;
+        if (oA < 0 || oA > nOther) continue;
+        const double prob = std::exp(lp + lcJ[r] + lcO[oA] - lt);
+        if (prob <= 0.0) continue;
+        const int s = nJ - r;
+        const double conc = c_p * choose2(s) + c_q * choose2(r);
+        aConc += prob * conc;
+        aDec += prob * (conc + double(p) * q * r * s);
+      }
+    }
+    eConc[M] = aConc;
+    eDec[M] = aDec;
+  }
+}
+
+// [[Rcpp::export]]
+List quartet_expect(const LogicalMatrix splits, const IntegerMatrix characters) {
+  const int n_splits = splits.ncol();
+  const int n_chars = characters.ncol();
+  NumericMatrix eConc(n_splits, n_chars), eDec(n_splits, n_chars);
+
+  for (int c = 0; c < n_chars; ++c) {
+    const CharCounts cc = char_counts(splits, characters, c);
+    const int nStates = int(cc.cnt.size());
+    if (nStates < 2) continue;
+    std::vector<double> ec, ed;
+    for (int a = 0; a < nStates - 1; ++a) {
+      for (int b = a + 1; b < nStates; ++b) {
+        expected_quartet_by_M(cc.cnt[a], cc.cnt[b], cc.tc, ec, ed);
+        for (int s = 0; s < n_splits; ++s) {
+          const int M = cc.mSideA[s];
+          eConc(s, c) += ec[M];
+          eDec(s, c) += ed[M];
+        }
+      }
+    }
+  }
+  return List::create(_["concordant"] = eConc, _["decisive"] = eDec);
+}
+
+// Exact E[m], E[m*A/wk], E[m*A/wc] for one NRQS state-pair over all side-A
+// sizes.  wc = (nI-1)+ (nJ-1)+ is fixed; wk = (mA-1)+ (tP-mA-1)+ and
+// m = min(wc, wk) vary with mA = p + r.  Mirrors R `.ExpectedTrit`.
+static void expected_nrqs_by_M(int nI, int nJ, int tc,
+                               std::vector<double>& eM,
+                               std::vector<double>& eMAwk,
+                               std::vector<double>& eMAwc) {
+  const int nOther = tc - nI - nJ;
+  const int tP = nI + nJ;
+  const double wc = posm1(nI) * posm1(nJ);
+  eM.assign(tc + 1, 0.0);
+  eMAwk.assign(tc + 1, 0.0);
+  eMAwc.assign(tc + 1, 0.0);
+  // Hoist the log-choose terms out of the M / p / r loops (the dominant cost).
+  std::vector<double> lcI(nI + 1), lcJ(nJ + 1), lcO(nOther + 1), lcT(tc + 1);
+  for (int p = 0; p <= nI; ++p) lcI[p] = R::lchoose(nI, p);
+  for (int r = 0; r <= nJ; ++r) lcJ[r] = R::lchoose(nJ, r);
+  for (int o = 0; o <= nOther; ++o) lcO[o] = R::lchoose(nOther, o);
+  for (int M = 0; M <= tc; ++M) lcT[M] = R::lchoose(tc, M);
+  for (int M = 0; M <= tc; ++M) {
+    const double lt = lcT[M];
+    double aM = 0.0, aWk = 0.0, aWc = 0.0;
+    const int pMax = std::min(nI, M);
+    for (int p = 0; p <= pMax; ++p) {
+      const int q = nI - p;
+      const double lp = lcI[p];
+      const int rMax = std::min(nJ, M - p);
+      for (int r = 0; r <= rMax; ++r) {
+        const int oA = M - p - r;
+        if (oA < 0 || oA > nOther) continue;
+        const double prob = std::exp(lp + lcJ[r] + lcO[oA] - lt);
+        if (prob <= 0.0) continue;
+        const int s = nJ - r;
+        const int mA = p + r;
+        const double A = posm1(p) * posm1(s) + posm1(q) * posm1(r);
+        const double wk = posm1(mA) * posm1(tP - mA);
+        const double m = std::min(wc, wk);
+        aM += prob * m;
+        if (wk > 0.0) aWk += prob * m * A / wk;
+        if (wc > 0.0) aWc += prob * m * A / wc;
+      }
+    }
+    eM[M] = aM;
+    eMAwk[M] = aWk;
+    eMAwc[M] = aWc;
+  }
+}
+
+// [[Rcpp::export]]
+List nrqs_expect(const LogicalMatrix splits, const IntegerMatrix characters) {
+  const int n_splits = splits.ncol();
+  const int n_chars = characters.ncol();
+  NumericMatrix numEdge(n_splits, n_chars);
+  NumericMatrix numChar(n_splits, n_chars);
+  NumericMatrix denM(n_splits, n_chars);
+
+  for (int c = 0; c < n_chars; ++c) {
+    const CharCounts cc = char_counts(splits, characters, c);
+    const int nStates = int(cc.cnt.size());
+    if (nStates < 2) continue;
+    std::vector<double> eM, eWk, eWc;
+    for (int a = 0; a < nStates - 1; ++a) {
+      for (int b = a + 1; b < nStates; ++b) {
+        expected_nrqs_by_M(cc.cnt[a], cc.cnt[b], cc.tc, eM, eWk, eWc);
+        for (int s = 0; s < n_splits; ++s) {
+          const int M = cc.mSideA[s];
+          denM(s, c) += eM[M];
+          numEdge(s, c) += eWk[M];
+          numChar(s, c) += eWc[M];
+        }
+      }
+    }
+  }
+  return List::create(_["numEdge"] = numEdge, _["numChar"] = numChar,
+                      _["denM"] = denM);
+}
diff --git a/src/expected_mi.cpp b/src/expected_mi.cpp
index c759a5547..289142f10 100644
--- a/src/expected_mi.cpp
+++ b/src/expected_mi.cpp
@@ -46,7 +46,7 @@ inline double l2factorial(int n) {
 //' block sizes (marginals) of each partition are fixed but the items are
 //' associated at random.  Subtracting this baseline from an observed mutual
 //' information yields a chance-corrected ("adjusted") mutual information, as
-//' applied by [`SiteConcordance`]`(normalize = TRUE)`; it is most material for
+//' applied by [`SiteConcordance`]`(chanceCorrect = TRUE)`; it is most material for
 //' small or unbalanced partitions, where raw mutual information is appreciably
 //' inflated by chance agreement.
 //'
diff --git a/tests/testthat/test-Concordance.R b/tests/testthat/test-Concordance.R
index 7d218dc4a..99a36bbdb 100644
--- a/tests/testthat/test-Concordance.R
+++ b/tests/testthat/test-Concordance.R
@@ -20,6 +20,11 @@ test_that("_Concordance() handles tip mismatch", {
 })
 
 test_that("QuartetConcordance() works", {
+  # This block locks the raw-quartet contract, so it names `unit` and
+  # `chanceCorrect` explicitly rather than relying on the (NRQS) defaults.
+  Naive <- function(...) {
+    QuartetConcordance(..., unit = "quartet", chanceCorrect = FALSE)
+  }
   tree <- BalancedTree(8)
   splits <- as.Splits(tree)
   mataset <- matrix(c(0, 0, 0, 0, 1, 1, 1, 1,  0,
@@ -32,15 +37,15 @@ test_that("QuartetConcordance() works", {
   expect_error(QuartetConcordance(tree, mataset),
                "`dataset` must be a phyDat object")
   dat <- MatrixToPhyDat(mataset)
-  expect_equal(unname(QuartetConcordance(tree, dat[, 1])), rep(1, 5))
+  expect_equal(unname(Naive(tree, dat[, 1])), rep(1, 5))
   # plot(tree); nodelabels();
-  expect_equal(QuartetConcordance(tree, dat[, 2]),
+  expect_equal(Naive(tree, dat[, 2]),
                c("10" = 1/9, "11" = 0, "12" = 0,
                  "13" = 1/9, "14" = 0, "15" = 0)[names(as.Splits(tree))])
   
   allQuartets <- combn(8, 4)
   for (charI in seq_len(ncol(mataset))) {
-    qc <- QuartetConcordance(tree, dat[, charI])
+    qc <- Naive(tree, dat[, charI])
     for (splitI in seq_along(splits)) {
       split <- splits[[splitI]]
       logiSplit <- as.logical(split)
@@ -68,7 +73,7 @@ test_that("QuartetConcordance() works", {
     }
   }
   
-  expect_equal(QuartetConcordance(tree, dat[, c(1:4, 6)]),
+  expect_equal(Naive(tree, dat[, c(1:4, 6)]),
                c("10" = (36 + 2 + 9 + 12) / (36 + 18 + 18 + 12 + 6),
                  "11" = ( 6 + 0 + 6 +  2) / ( 6 +  9 +  6 +  2 + 1),
                  "12" = ( 6 + 0 + 0 +  2) / ( 6 +  9 +  9 +  2 + 1),
@@ -133,6 +138,267 @@ test_that("QuartetConcordance() handles non-integer data", {
                QuartetConcordance(tree, MatrixToPhyDat(intSet)))
 })
 
+test_that("QuartetConcordance() unit = 'nrqs' locks the contract", {
+  tree <- BalancedTree(8)
+
+  # Character identical to the {t1..t4 | t5..t8} split: that split scores 1,
+  # and it is the *only* split scoring 1 (nested/crossing get partial credit).
+  identChar <- MatrixToPhyDat(matrix(
+    c(0, 0, 0, 0, 1, 1, 1, 1), 8, dimnames = list(paste0("t", 1:8), NULL)))
+  qt <- QuartetConcordance(tree, identChar, unit = "nrqs",
+                           chanceCorrect = FALSE)
+  expect_equal(max(qt, na.rm = TRUE), 1)
+  expect_equal(sum(abs(qt - 1) < 1e-9, na.rm = TRUE), 1L)
+  expect_true(all(qt >= 0 & qt <= 1, na.rm = TRUE))
+
+  # `unit = "nrqs"` and `chanceCorrect = TRUE` are the defaults, so a bare call
+  # is byte-identical to naming them.
+  expect_identical(QuartetConcordance(tree, identChar),
+                   QuartetConcordance(tree, identChar, unit = "nrqs",
+                                      chanceCorrect = TRUE))
+  # Coverage normalisation makes NRQS no laxer than quartet.
+  expect_true(mean(qt, na.rm = TRUE) <=
+                mean(QuartetConcordance(tree, identChar, unit = "quartet",
+                                        chanceCorrect = FALSE),
+                     na.rm = TRUE))
+
+  # Uninformative characters (constant / autapomorphy) carry no NRQS -> NA.
+  autap <- MatrixToPhyDat(matrix(
+    c(0, 0, 0, 0, 0, 0, 0, 1), 8, dimnames = list(paste0("t", 1:8), NULL)))
+  expect_equal(unname(QuartetConcordance(tree, autap, unit = "nrqs")),
+               rep(NA_real_, 5))
+
+  # `unit` is validated.
+  expect_error(QuartetConcordance(tree, identChar, unit = "trits"),
+               "should be one of")
+
+  # `return` is matched partially, as on the quartet path.
+  ByReturn <- function(r) {
+    QuartetConcordance(tree, identChar, return = r, unit = "nrqs")
+  }
+  expect_equal(ByReturn("char"), ByReturn("ch"))
+  expect_error(ByReturn("site"), "must .* match")
+})
+
+test_that("QuartetConcordance() unit = 'nrqs' gives nested partial credit", {
+  # {t1..t5 | t6,t7,t8} split; state {t1,t2,t4} nested within side A but not a
+  # clade, giving cells (p,q,r,s) = (3,0,2,3) -> A/Wk = 4/8 = 0.5 exactly.
+  tree <- ape::read.tree(text = "(((((t1,t2),t3),t4),t5),(t6,(t7,t8)));")
+  char <- MatrixToPhyDat(matrix(
+    c(0, 0, 1, 0, 1, 1, 1, 1), 8, dimnames = list(paste0("t", 1:8), NULL)))
+  qt <- QuartetConcordance(tree, char, unit = "nrqs", chanceCorrect = FALSE)
+
+  # The split isolating {t6,t7,t8} scores exactly the nested value 0.5.
+  sp <- as.Splits(tree)
+  tips <- TipLabels(sp)
+  member <- vapply(seq_along(sp), function(i) {
+    side <- tips[as.logical(sp[[i]])]
+    setequal(side, c("t6", "t7", "t8")) ||
+      setequal(setdiff(tips, side), c("t6", "t7", "t8"))
+  }, logical(1))
+  expect_equal(unname(qt[member]), 0.5)
+})
+
+test_that("QuartetConcordance() unit = 'nrqs' supports multistate", {
+  # Multistate no longer errors: NRQS sum over state-pairs (same currency).
+  tree <- BalancedTree(8)
+  ms <- MatrixToPhyDat(matrix(
+    c(0, 0, 1, 1, 2, 2, 2, 0,      # 3 states
+      0, 0, 0, 1, 1, 2, 3, 3), 8,  # 4 states
+    dimnames = list(paste0("t", 1:8), NULL)))
+  qt <- QuartetConcordance(tree, ms, unit = "nrqs", chanceCorrect = FALSE)
+  expect_length(qt, 5L)
+  expect_true(all(qt >= 0 & qt <= 1, na.rm = TRUE))
+  expect_false(anyNA(qt))  # both characters are informative on this tree
+})
+
+test_that("QuartetConcordance() unit = 'nrqs' scores 1 iff split displayed", {
+  # A multistate character need not be *identical* to a split to score 1: it
+  # scores full marks for every split its own tree displays (each state block
+  # wholly on one side; the split-orthogonal state-pair drops via M = 0).  This
+  # generalises "only identical scores 1" from binary to multistate -- a cleanly
+  # refining reproductive character fully supports the clade it refines.
+  tree <- ape::read.tree(text = "((((t1,t2),t3),((t4,t5),t6)),((t7,t8),t9));")
+  char <- MatrixToPhyDat(matrix(
+    c(0, 0, 0, 1, 1, 1, 2, 2, 2), 9,
+    dimnames = list(paste0("t", 1:9), NULL)))
+  qt <- QuartetConcordance(tree, char, unit = "nrqs", chanceCorrect = FALSE)
+  sp <- as.Splits(tree)
+  tips <- TipLabels(sp)
+  atSplit <- function(members) {
+    i <- which(vapply(seq_along(sp), function(j) {
+      side <- tips[as.logical(sp[[j]])]
+      setequal(side, members) || setequal(setdiff(tips, side), members)
+    }, logical(1)))
+    unname(qt[i])
+  }
+  # Splits the character displays -> exactly 1
+  expect_equal(atSplit(c("t1", "t2", "t3")), 1)
+  expect_equal(atSplit(c("t4", "t5", "t6")), 1)
+  expect_equal(atSplit(c("t1", "t2", "t3", "t4", "t5", "t6")), 1)
+  # A split cutting through a state block is only partially supported
+  expect_true(atSplit(c("t1", "t2")) < 1)
+})
+
+test_that("QuartetConcordance() 'nrqs' return = 'char' has no ceiling of 1", {
+  # Documented limitation, not an oversight: the character score is the
+  # M-weighted mean across EVERY split, and a character agrees exactly with at
+  # most one of them while being merely compatible with -- silent about -- the
+  # rest.  Coverage scores silence as failure to cover, so even a character
+  # identical to a split of the tree falls short of 1.  `unit = "quartet"` is
+  # the per-character measure that does attain 1.
+  tree <- BalancedTree(8)
+  tips <- TipLabels(tree)
+  sp <- as.Splits(tree)
+  ByUnit <- function(u, r) {
+    QuartetConcordance(tree, char, return = r, unit = u, chanceCorrect = FALSE)
+  }
+  for (k in seq_along(sp)) {
+    char <- MatrixToPhyDat(matrix(
+      as.integer(as.logical(sp[[k]])), length(tips),
+      dimnames = list(tips, NULL)))
+
+    # Per-split bound holds: the character scores 1 at the split it duplicates
+    expect_equal(max(ByUnit("nrqs", "edge")), 1)
+
+    # Naive currency reaches 1 on the character path ...
+    expect_equal(ByUnit("quartet", "char")[[1]], 1)
+
+    # ... the NRQS currency does not, and must not be read against 1
+    nrqsChar <- ByUnit("nrqs", "char")[[1]]
+    expect_gt(nrqsChar, 0)
+    expect_lt(nrqsChar, 1)
+  }
+})
+
+test_that("QuartetConcordance() unit = 'nrqs' chanceCorrect contract", {
+  tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));")
+  ident <- MatrixToPhyDat(matrix(
+    c(0, 0, 0, 0, 0, 0, 1, 1, 1), 9,
+    dimnames = list(paste0("t", 1:9), NULL)))
+
+  # chanceCorrect = TRUE is the default; re-zeroing can only lower the measure,
+  # since a non-negative baseline is subtracted.
+  expect_identical(
+    QuartetConcordance(tree, ident, unit = "nrqs", chanceCorrect = TRUE),
+    QuartetConcordance(tree, ident, unit = "nrqs"))
+  expect_true(all(
+    QuartetConcordance(tree, ident, unit = "nrqs", chanceCorrect = TRUE) <=
+      QuartetConcordance(tree, ident, unit = "nrqs", chanceCorrect = FALSE),
+    na.rm = TRUE))
+
+  # Invalid chanceCorrect is rejected.
+  for (bad in list(0, -3, "x")) {
+    expect_error(QuartetConcordance(tree, ident, chanceCorrect = bad),
+                 "positive integer")
+  }
+
+  # Chance correction also works for the raw quartet currency.
+  qc <- QuartetConcordance(tree, ident, unit = "quartet", chanceCorrect = TRUE)
+  expect_type(qc, "double")
+  expect_false(anyNA(qc))
+
+  # An identical (displayed) split still scores exactly 1 after re-zeroing:
+  # .Rezero(1, z) == 1, so chance correction lifts the floor without moving the
+  # ceiling.
+  sp <- as.Splits(tree)
+  tips <- TipLabels(sp)
+  col789 <- which(vapply(seq_along(sp), function(j) {
+    side <- tips[as.logical(sp[[j]])]
+    setequal(side, c("t7", "t8", "t9")) ||
+      setequal(setdiff(tips, side), c("t7", "t8", "t9"))
+  }, logical(1)))
+  corrected <- QuartetConcordance(tree, ident, unit = "nrqs",
+                                  chanceCorrect = TRUE)
+  expect_equal(unname(corrected[col789]), 1)
+})
+
+test_that("QuartetConcordance() unit = 'nrqs' exact baseline matches MC", {
+  tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));")
+  dat <- MatrixToPhyDat(matrix(
+    c(0, 0, 0, 0, 0, 0, 1, 1, 1,      # displays {t7,t8,t9}
+      0, 0, 0, 0, 0, 0, 0, 1, 1,      # nested: {t8,t9}
+      0, 0, 1, 1, 0, 0, 1, 1, 0,      # crossing (binary)
+      0, 0, 0, 1, 1, 1, 2, 2, 2), 9,  # 3-state: exercises the random-wk path
+    dimnames = list(paste0("t", 1:9), NULL)))
+
+  # The exact hypergeometric baseline and the Monte-Carlo tip-shuffle baseline
+  # target the same null, so they agree to within MC error -- for edge and char,
+  # under both weightings, and crucially through the MULTISTATE pair path where a
+  # pair's wk and M = min(w_c, w_k) are themselves random under the null.
+  for (ret in c("edge", "char")) {
+    for (w in c(TRUE, FALSE)) {
+      exact <- QuartetConcordance(tree, dat, unit = "nrqs", return = ret,
+                                  weight = w, chanceCorrect = TRUE)
+      set.seed(1)
+      mc <- QuartetConcordance(tree, dat, unit = "nrqs", return = ret,
+                               weight = w, chanceCorrect = 8000)
+      # Same cells resolve to NA under exact and MC (guards the cell-matching in
+      # the weight = FALSE path), and the finite values agree within MC error.
+      expect_identical(is.na(exact), is.na(mc))
+      expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04)
+    }
+  }
+})
+
+test_that("QuartetConcordance() unit = 'quartet' exact baseline matches MC", {
+  tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));")
+  dat <- MatrixToPhyDat(matrix(
+    c(0, 0, 0, 0, 0, 0, 1, 1, 1,      # displays {t7,t8,t9}
+      0, 0, 0, 0, 0, 0, 0, 1, 1,      # nested
+      0, 0, 1, 1, 0, 0, 1, 1, 0,      # crossing (binary)
+      0, 0, 0, 1, 1, 1, 2, 2, 2), 9,  # 3-state: multistate decisive path
+    dimnames = list(paste0("t", 1:9), NULL)))
+
+  # The raw-currency chance baseline (exact E[conc]/E[dec] from the
+  # hypergeometric pmf) matches the Monte-Carlo tip-shuffle baseline within MC
+  # error, across edge/char and both weightings, incl. the multistate path.
+  for (ret in c("edge", "char")) {
+    for (w in c(TRUE, FALSE)) {
+      exact <- QuartetConcordance(tree, dat, unit = "quartet", return = ret,
+                                  weight = w, chanceCorrect = TRUE)
+      set.seed(1)
+      mc <- QuartetConcordance(tree, dat, unit = "quartet", return = ret,
+                               weight = w, chanceCorrect = 8000)
+      expect_identical(is.na(exact), is.na(mc))
+      expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04)
+    }
+  }
+
+  # The published raw measure is recovered by naming both non-default settings,
+  # and the default no longer returns it.
+  for (ret in c("edge", "char")) for (w in c(TRUE, FALSE)) {
+    raw <- QuartetConcordance(tree, dat, unit = "quartet", return = ret,
+                              weight = w, chanceCorrect = FALSE)
+    expect_true(all(raw >= 0 & raw <= 1, na.rm = TRUE))
+    expect_false(identical(
+      raw,
+      QuartetConcordance(tree, dat, return = ret, weight = w)))
+  }
+})
+
+test_that("QuartetConcordance() 'nrqs' correction: conflict below random", {
+  tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));")
+  support <- MatrixToPhyDat(matrix(   # agrees with the tree: displays {t7,t8,t9}
+    c(0, 0, 0, 0, 0, 0, 1, 1, 1), 9,
+    dimnames = list(paste0("t", 1:9), NULL)))
+  conflict <- MatrixToPhyDat(matrix(  # crosses the two halves of the tree
+    c(0, 0, 1, 1, 0, 0, 1, 1, 0), 9,
+    dimnames = list(paste0("t", 1:9), NULL)))
+
+  sChar <- QuartetConcordance(tree, support, return = "char", unit = "nrqs",
+                              chanceCorrect = TRUE)
+  cNone <- QuartetConcordance(tree, conflict, return = "char", unit = "nrqs",
+                              chanceCorrect = FALSE)
+  cChar <- QuartetConcordance(tree, conflict, return = "char", unit = "nrqs",
+                              chanceCorrect = TRUE)
+  # A character that agrees with the tree scores above random expectation;
+  # a crossing character is pulled below it, and below its uncorrected score.
+  expect_gt(sChar, 0)
+  expect_lt(cChar, cNone)
+  expect_lt(cChar, 0)
+})
+
 test_that(".Rezero() works", {
   expect_equal(TreeSearch:::.Rezero(seq(0, 1, by = 0.1), 0.1), -1:9 / 9)
 })
@@ -218,13 +484,15 @@ test_that("ClusteringConcordance() gives sensible values", {
   randomset <- matrix(sample(0:1, 8 * 1000, replace = TRUE), 8,
                       dimnames = list(letters[1:8], NULL))
   rat <- MatrixToPhyDat(randomset)
-  expect_equal(ClusteringConcordance(tree, rat, normalize = TRUE), c("10" = 0),
+  expect_equal(ClusteringConcordance(tree, rat, chanceCorrect = TRUE),
+               c("10" = 0),
                tolerance = 0.05)
 })
 
 test_that("ClusteringConcordance(return = 'char') Monte-Carlo handles ambiguity", {
   # Regression for T-330: characters whose ambiguous tokens drop different tips
-  # give `charSplits` over heterogeneous tip sets. The Monte-Carlo `normalize`
+  # give `charSplits` over heterogeneous tip sets. The Monte-Carlo
+  # `chanceCorrect`
   # path scored a *list* of random trees against that list in one call, which
   # could not reconcile a common label set ("Old and new labels must match").
   tree <- ape::read.tree(text = "((a, b, c, d, e), (f, g, h));")
@@ -238,7 +506,7 @@ test_that("ClusteringConcordance(return = 'char') Monte-Carlo handles ambiguity"
 
   set.seed(1)
   # Previously errored: "Old and new labels must match"
-  cc <- ClusteringConcordance(tree, dat, return = "char", normalize = 10L)
+  cc <- ClusteringConcordance(tree, dat, return = "char", chanceCorrect = 10L)
   nChar <- length(attr(dat, "index"))
   expect_length(cc, nChar)
   expect_length(attr(cc, "mcse"), nChar)
@@ -246,7 +514,8 @@ test_that("ClusteringConcordance(return = 'char') Monte-Carlo handles ambiguity"
   expect_equal(unname(cc[2]), 1)
   expect_equal(unname(attr(cc, "mcse")[2]), 0)
   # Monte-Carlo score never exceeds the un-normalized score (subtracts a baseline)
-  bare <- ClusteringConcordance(tree, dat, return = "char", normalize = FALSE)
+  bare <- ClusteringConcordance(tree, dat, return = "char",
+                                chanceCorrect = FALSE)
   finite <- is.finite(cc) & is.finite(bare)
   expect_true(all(cc[finite] <= bare[finite] + 1e-8))
 })
@@ -387,7 +656,7 @@ test_that("QuartetConcordance() handles a {0,-} contrast level (#108)", {
   dat <- MatrixToPhyDat(m)
   expect_no_error(ret <- QuartetConcordance(tree, dat))
   expect_length(ret, length(as.Splits(tree)))
-  expect_true(all(is.na(ret) | (ret >= 0 & ret <= 1)))
+  expect_true(all(is.na(ret) | (is.finite(ret) & ret <= 1)))
 
   # `{0-}` must convey no grouping information, exactly like `?` -- not
   # merely avoid crashing.
diff --git a/vignettes/tree-search.Rmd b/vignettes/tree-search.Rmd
index cf0694bcd..1cd83b07c 100644
--- a/vignettes/tree-search.Rmd
+++ b/vignettes/tree-search.Rmd
@@ -175,6 +175,10 @@ smaller number of characters contradicting it.
 whether a dataset is unanimous or divided in support of a grouping, 
 independently of the method of tree reconstruction.
 
+By default, these are corrected for chance: 1 denotes unanimous support, 0 no
+more agreement than would be expected if the character's tokens were shuffled
+across the leaves, and negative values conflict.
+
 ```{r concordance}
 concordance <- QuartetConcordance(strict, vinther)
 
@@ -185,7 +189,7 @@ concordance <- QuartetConcordance(strict, vinther)
 par(mar = rep(0, 4), cex = 0.8)
 plot(strict)
 TreeTools::LabelSplits(strict, signif(concordance, 3),
-                       col = TreeTools::SupportColor(concordance / max(concordance)),
+                       col = TreeTools::SupportColor(pmax(0, concordance)),
                        frame = "none", pos = 3L)
 ```