diff --git a/R/Concordance.R b/R/Concordance.R
index 2ba797245..51a484d22 100644
--- a/R/Concordance.R
+++ b/R/Concordance.R
@@ -87,15 +87,29 @@ NULL
 #' 
 #'   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 normalize Controls the zero point of the concordance scale by
+#' subtracting the value expected under a chance (fixed-marginal) null, in which
+#' each character's tokens are reassigned at random across the leaves while its
+#' state frequencies and the split sizes are held fixed.
+#'   - `FALSE`: no chance correction; the measure is scaled only by its maximum,
+#'     so 1 marks a perfect match and 0 the measure's own (typically positive)
+#'     floor.
+#'   - `TRUE`: subtract the expected value under the null, so that 0 marks random
+#'     expectation and negative values fall below it. `ClusteringConcordance()`
+#'     uses an analytical approximation to the expected mutual information (fast
+#'     and generally accurate for large trees, ~200+ taxa, but neglecting
+#'     correlation between splits); `QuartetConcordance()` uses the exact
+#'     hypergeometric expectation (for either `unit`).
+#'   - a positive integer `n`: estimate that expectation by Monte Carlo.
+#'     `ClusteringConcordance()` fits each character to `n` random trees (and
+#'     returns Monte-Carlo standard errors, more accurate for small trees where
+#'     the analytical approximation is biased); `QuartetConcordance()` averages
+#'     over `n` random reassignments of each character's tokens.
+#'
+#'   In all cases 1 corresponds to the maximum attainable value.  For
+#'   `QuartetConcordance()` (either `unit`), chance-corrected values are returned
+#'   unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before
+#'   plotting with [QCol()] / [QACol()].
 #' 
 #' @returns
 #' `ClusteringConcordance(return = "all")` returns a 3D array where each
@@ -765,8 +779,31 @@ MutualClusteringConcordance <- function(tree, dataset) {
 #' `QuartetConcordance(return = "char")` returns a numeric vector giving the
 #' concordance index calculated at each site, averaged across all splits.
 #'
+#' With `unit = "trit"` (see below) the same vectors are returned, but scored in
+#' redundancy-corrected currency: values are typically lower, and reach 1 only
+#' where a split is *displayed* by the character rather than merely concordant
+#' with many of its quartets.
+#'
 #' @param weight Logical specifying whether to weight sites according to the
 #' number of quartets they are decisive for.
+#' @param unit Character specifying the currency in which quartets are counted:
+#'   - `"quartet"` (default): each resolved quartet counts once, so a character
+#'     is credited with the full combinatorial volume of quartets it resolves;
+#'   - `"trit"`: quartets are counted as the *independent* information they carry
+#'     \insertCite{Nelson1992}{TreeSearch}.  Because the quartets resolved by a
+#'     split of sizes \eqn{(k, t - k)} are logically redundant --
+#'     \eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} -- only
+#'     \eqn{(k - 1)(t - k - 1)} of the \eqn{\binom{k}{2}\binom{t - k}{2}}
+#'     resolved quartets are independent.  Concordance is then measured against
+#'     this reduced (redundancy-corrected) content, so that only a character
+#'     whose split is *identical* to the tree split scores full marks; nested
+#'     (compatible) characters receive genuine partial support, and crossing
+#'     (incompatible) characters score lower still.  A multistate character is
+#'     scored in the same currency: trits are counted independently within each
+#'     pair of states (sizes \eqn{n_i}, \eqn{n_j} give a per-state-pair weight of
+#'     \eqn{4 / (n_i n_j)}) and summed across pairs, so multistate and binary
+#'     characters remain directly comparable.
+#' @references \insertAllCited{}
 #' @importFrom ape keep.tip
 #' @importFrom cli cli_progress_bar cli_progress_update
 #' @importFrom utils combn
@@ -776,7 +813,9 @@ QuartetConcordance <- function(
   tree,
   dataset = NULL,
   weight = TRUE,
-  return = "edge"
+  return = "edge",
+  unit = c("quartet", "trit"),
+  normalize = FALSE
 ) {
   if (is.null(dataset)) {
     warning("Cannot calculate concordance without `dataset`.")
@@ -785,6 +824,15 @@ QuartetConcordance <- function(
   if (!inherits(dataset, "phyDat")) {
     stop("`dataset` must be a phyDat object.")
   }
+  unit <- match.arg(unit)
+  if (!isFALSE(normalize)) {
+    if (!isTRUE(normalize)) {
+      if (!is.numeric(normalize) || length(normalize) != 1L ||
+          is.na(normalize) || normalize < 1) {
+        stop("`normalize` 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.")
@@ -816,14 +864,37 @@ QuartetConcordance <- function(
     dimnames = dimnames(characters)
   )
   
+  options <- c("character", "site", "default")
+  return <- options[[pmatch(tolower(trimws(return)), options,
+                            nomatch = length(options))]]
+
+  if (unit == "trit") {
+    # Return:
+    return(.TritConcordance(logiSplits, charInt, weight, return, splits,
+                            normalize))
+  }
+
   raw_counts <- quartet_concordance(logiSplits, charInt)
 
   num <- raw_counts$concordant
   den <- raw_counts$decisive
-  options <- c("character", "site", "default")
-  return <- options[[pmatch(tolower(trimws(return)), options,
-                            nomatch = length(options))]]
-  
+
+  # Chance correction (option): re-zero the observed concordant/decisive ratio
+  # against the ratio expected under the same fixed-marginal null used for the
+  # trit 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 trits)
+  # or by Monte-Carlo tip-shuffle, then re-zero the pooled ratio.
+  doNorm <- !isFALSE(normalize)
+  if (doNorm) {
+    base <- if (isTRUE(normalize)) {
+      .QuartetExpect(charInt, logiSplits)
+    } else {
+      .QuartetMC(charInt, logiSplits, normalize)
+    }
+    eNum <- base[["concordant"]]
+    eDen <- base[["decisive"]]
+  }
 
   if (return == "default") {
     if (isTRUE(weight)) {
@@ -836,39 +907,334 @@ 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"
-    p <- num / den
     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 ("trit") 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 independent ("trits").
+#
+# A multistate 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 trits 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 trits 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 scores 1.  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.
+.TritConcordance <- function(logiSplits, charInt, weight, return, splits,
+                             normalize = FALSE) {
+  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 trit content (split-free)
+  doNorm <- !isFALSE(normalize)
+  if (doNorm) {
+    baseNumEdge <- baseNumChar <- baseDenM <- matrix(0, nSplit, nChar)
+    if (isTRUE(normalize)) {
+      # Exact hypergeometric expectation for all characters at once (C++;
+      # `trit_expect` in src/concordance_expect.cpp).  Uninformative characters
+      # yield zero pools, matching the observed loop's `wcTot > 0` guard.
+      be <- trit_expect(logiSplits, charInt)
+      baseNumEdge <- be[["numEdge"]]
+      baseNumChar <- be[["numChar"]]
+      baseDenM <- be[["denM"]]
+    }
+  }
+
+  for (ci in seq_len(nChar)) {
+    col <- charInt[, ci]
+    obs <- .CharTritContrib(col, logiSplits, pos)
+    numEdge[, ci] <- obs[["numEdge"]]
+    numChar[, ci] <- obs[["numChar"]]
+    denM[, ci] <- obs[["denM"]]
+    wcTot[ci] <- obs[["wcTot"]]
+    if (doNorm && !isTRUE(normalize) && obs[["wcTot"]] > 0) {
+      base <- .CharTritMC(col, logiSplits, pos, normalize)  # Monte-Carlo shuffle
+      baseNumEdge[, ci] <- base[["numEdge"]]
+      baseNumChar[, ci] <- base[["numChar"]]
+      baseDenM[, ci] <- base[["denM"]]
+    }
+  }
+
+  informative <- wcTot > 0
+
+  if (return == "default") {
+    # 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 trits 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 multistate 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 trit 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.
+.CharTritContrib <- 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 trits; 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 trit expectation E[m], E[m*A/wk], E[m*A/wc] is computed in C++
+# (`trit_expect` in src/concordance_expect.cpp) for all characters at once, and
+# consumed directly by `.TritConcordance`; there is no per-character R exact
+# helper.  The Monte-Carlo baseline below is the opt-in `normalize = <int>` path.
+
+# Monte-Carlo baseline: average `.CharTritContrib` 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.
+.CharTritMC <- 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 <- .CharTritContrib(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())
 
 # @param a must be a vector of length <= 2
diff --git a/R/RcppExports.R b/R/RcppExports.R
index afce48989..4d29eae1d 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)
+}
+
+trit_expect <- function(splits, characters) {
+    .Call(`_TreeSearch_trit_expect`, splits, characters)
+}
+
 #' Expected mutual information between two partitions
 #'
 #' Computes the mutual information expected purely by chance between two
diff --git a/dev/benchmarks/frac-quart/chance_baseline.R b/dev/benchmarks/frac-quart/chance_baseline.R
new file mode 100644
index 000000000..93b9d9690
--- /dev/null
+++ b/dev/benchmarks/frac-quart/chance_baseline.R
@@ -0,0 +1,86 @@
+# Validate the exact hypergeometric chance baseline used by
+# QuartetConcordance(unit = "trit", normalize = 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..0e8afba7f
--- /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
+# `trit_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 <- trit_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..b0f9142e9
--- /dev/null
+++ b/dev/benchmarks/frac-quart/validate_trit.R
@@ -0,0 +1,229 @@
+# Validation of the unit = "trit" 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 = "trit")
+  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 = "trit")
+  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 = "trit")
+  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 = "trit")
+e2 <- QuartetConcordance(refTree, binDat, return = "default",   unit = "trit")
+c1 <- QuartetConcordance(refTree, binDat, return = "char",      unit = "trit")
+c2 <- QuartetConcordance(refTree, binDat, return = "character", unit = "trit")
+c3 <- QuartetConcordance(refTree, binDat, return = "site",      unit = "trit")
+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")
+tEdge <- QuartetConcordance(refTree, binDat, unit = "trit")
+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/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd
index 9d6a0831b..b420da41c 100644
--- a/man/ConcordanceTable.Rd
+++ b/man/ConcordanceTable.Rd
@@ -36,17 +36,31 @@ 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{normalize}{Controls the zero point of the concordance scale by
+subtracting the value expected under a chance (fixed-marginal) null, in which
+each character's tokens are reassigned at random across the leaves while its
+state frequencies and the split sizes are held fixed.
 \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.
+\item \code{FALSE}: no chance correction; the measure is scaled only by its maximum,
+so 1 marks a perfect match and 0 the measure's own (typically positive)
+floor.
+\item \code{TRUE}: subtract the expected value under the null, so that 0 marks random
+expectation and negative values fall below it. \code{ClusteringConcordance()}
+uses an analytical approximation to the expected mutual information (fast
+and generally accurate for large trees, ~200+ taxa, but neglecting
+correlation between splits); \code{QuartetConcordance()} uses the exact
+hypergeometric expectation (for either \code{unit}).
+\item a positive integer \code{n}: estimate that expectation by Monte Carlo.
+\code{ClusteringConcordance()} fits each character to \code{n} random trees (and
+returns Monte-Carlo standard errors, more accurate for small trees where
+the analytical approximation is biased); \code{QuartetConcordance()} averages
+over \code{n} random reassignments of each character's tokens.
 }
 
-In all cases, 1 corresponds to the maximal attainable MI for the pair
-(\code{hBest}), and 0 corresponds to the chosen expectation.}
+In all cases 1 corresponds to the maximum attainable value.  For
+\code{QuartetConcordance()} (either \code{unit}), chance-corrected values are returned
+unclamped (they may fall below \eqn{-1}); clamp 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 7dc208b54..f8b3c6cc0 100644
--- a/man/SiteConcordance.Rd
+++ b/man/SiteConcordance.Rd
@@ -13,7 +13,14 @@ ClusteringConcordance(tree, dataset, return = "edge", normalize = TRUE)
 
 MutualClusteringConcordance(tree, dataset)
 
-QuartetConcordance(tree, dataset = NULL, weight = TRUE, return = "edge")
+QuartetConcordance(
+  tree,
+  dataset = NULL,
+  weight = TRUE,
+  return = "edge",
+  unit = c("quartet", "trit"),
+  normalize = FALSE
+)
 
 PhylogeneticConcordance(tree, dataset)
 
@@ -37,20 +44,54 @@ 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{normalize}{Controls the zero point of the concordance scale by
+subtracting the value expected under a chance (fixed-marginal) null, in which
+each character's tokens are reassigned at random across the leaves while its
+state frequencies and the split sizes are held fixed.
 \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.
+\item \code{FALSE}: no chance correction; the measure is scaled only by its maximum,
+so 1 marks a perfect match and 0 the measure's own (typically positive)
+floor.
+\item \code{TRUE}: subtract the expected value under the null, so that 0 marks random
+expectation and negative values fall below it. \code{ClusteringConcordance()}
+uses an analytical approximation to the expected mutual information (fast
+and generally accurate for large trees, ~200+ taxa, but neglecting
+correlation between splits); \code{QuartetConcordance()} uses the exact
+hypergeometric expectation (for either \code{unit}).
+\item a positive integer \code{n}: estimate that expectation by Monte Carlo.
+\code{ClusteringConcordance()} fits each character to \code{n} random trees (and
+returns Monte-Carlo standard errors, more accurate for small trees where
+the analytical approximation is biased); \code{QuartetConcordance()} averages
+over \code{n} random reassignments of each character's tokens.
 }
 
-In all cases, 1 corresponds to the maximal attainable MI for the pair
-(\code{hBest}), and 0 corresponds to the chosen expectation.}
+In all cases 1 corresponds to the maximum attainable value.  For
+\code{QuartetConcordance()} (either \code{unit}), chance-corrected values are returned
+unclamped (they may fall below \eqn{-1}); clamp 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.}
+
+\item{unit}{Character specifying the currency in which quartets are counted:
+\itemize{
+\item \code{"quartet"} (default): each resolved quartet counts once, so a character
+is credited with the full combinatorial volume of quartets it resolves;
+\item \code{"trit"}: quartets are counted as the \emph{independent} information they carry
+\insertCite{Nelson1992}{TreeSearch}.  Because the quartets resolved by a
+split of sizes \eqn{(k, t - k)} are logically redundant --
+\eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} -- only
+\eqn{(k - 1)(t - k - 1)} of the \eqn{\binom{k}{2}\binom{t - k}{2}}
+resolved quartets are independent.  Concordance is then measured against
+this reduced (redundancy-corrected) content, so that only a character
+whose split is \emph{identical} to the tree split scores full marks; nested
+(compatible) characters receive genuine partial support, and crossing
+(incompatible) characters score lower still.  A multistate character is
+scored in the same currency: trits are counted independently within each
+pair of states (sizes \eqn{n_i}, \eqn{n_j} give a per-state-pair weight of
+\eqn{4 / (n_i n_j)}) and summed across pairs, so multistate and binary
+characters remain directly comparable.
+}}
 }
 \value{
 \code{ClusteringConcordance(return = "all")} returns a 3D array where each
@@ -99,6 +140,11 @@ each corresponding split in \code{tree}.
 \code{QuartetConcordance(return = "char")} returns a numeric vector giving the
 concordance index calculated at each site, averaged across all splits.
 
+With \code{unit = "trit"} (see below) the same vectors are returned, but scored in
+redundancy-corrected currency: values are typically lower, and reach 1 only
+where a split is \emph{displayed} by the character rather than merely concordant
+with many of its quartets.
+
 \code{PhylogeneticConcordance()} returns a numeric vector giving the
 phylogenetic information of each split in \code{tree}, named according to the
 split's internal numbering.
@@ -236,6 +282,8 @@ myMatrix <- congreveLamsdellMatrices[[10]]
 ClusteringConcordance(TreeTools::NJTree(myMatrix), myMatrix)
 }
 \references{
+\insertAllCited{}
+
 \insertAllCited{}
 }
 \seealso{
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index dc52e9b59..f46299b2e 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
+}
+// trit_expect
+List trit_expect(const LogicalMatrix splits, const IntegerMatrix characters);
+RcppExport SEXP _TreeSearch_trit_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(trit_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 7b166a616..890fd4255 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_trit_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);
@@ -75,6 +77,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_trit_expect",       (DL_FUNC) &_TreeSearch_trit_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..4faefac8e
--- /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(normalize = 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 trit 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_trit_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 trit_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_trit_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/tests/testthat/test-Concordance.R b/tests/testthat/test-Concordance.R
index 05a75f53f..b455f6012 100644
--- a/tests/testthat/test-Concordance.R
+++ b/tests/testthat/test-Concordance.R
@@ -133,6 +133,226 @@ test_that("QuartetConcordance() handles non-integer data", {
                QuartetConcordance(tree, MatrixToPhyDat(intSet)))
 })
 
+test_that("QuartetConcordance() unit = 'trit' 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 = "trit")
+  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 = "quartet" is the default and is byte-identical to omitting `unit`.
+  expect_identical(QuartetConcordance(tree, identChar),
+                   QuartetConcordance(tree, identChar, unit = "quartet"))
+  # Coverage normalisation makes trit no laxer than quartet.
+  expect_true(mean(qt, na.rm = TRUE) <=
+                mean(QuartetConcordance(tree, identChar, unit = "quartet"),
+                     na.rm = TRUE))
+
+  # Uninformative characters (constant / autapomorphy) carry no trits -> 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 = "trit")),
+               rep(NA_real_, 5))
+
+  # `unit` is validated.
+  expect_error(QuartetConcordance(tree, identChar, unit = "trits"),
+               "should be one of")
+
+  # `return` aliases mirror the quartet path.
+  expect_equal(QuartetConcordance(tree, identChar, return = "edge", unit = "trit"),
+               QuartetConcordance(tree, identChar, return = "default", unit = "trit"))
+  cA <- QuartetConcordance(tree, identChar, return = "char", unit = "trit")
+  expect_equal(cA,
+               QuartetConcordance(tree, identChar, return = "character", unit = "trit"))
+  expect_equal(cA,
+               QuartetConcordance(tree, identChar, return = "site", unit = "trit"))
+})
+
+test_that("QuartetConcordance() unit = 'trit' 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 = "trit")
+
+  # 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 = 'trit' supports multistate", {
+  # Multistate no longer errors: trits 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 = "trit")
+  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 = 'trit' 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 = "trit")
+  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() unit = 'trit' normalize 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)))
+
+  # normalize = FALSE is the default and leaves the measure untouched.
+  expect_identical(
+    QuartetConcordance(tree, ident, unit = "trit", normalize = FALSE),
+    QuartetConcordance(tree, ident, unit = "trit"))
+
+  # Invalid normalize is rejected.
+  expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = 0),
+               "positive integer")
+  expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = -3),
+               "positive integer")
+  expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = "x"),
+               "positive integer")
+
+  # Chance correction also works for the raw quartet currency.
+  qc <- QuartetConcordance(tree, ident, unit = "quartet", normalize = 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 = "trit", normalize = TRUE)
+  expect_equal(unname(corrected[col789]), 1)
+})
+
+test_that("QuartetConcordance() unit = 'trit' 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 = "trit", return = ret,
+                                  weight = w, normalize = TRUE)
+      set.seed(1)
+      mc <- QuartetConcordance(tree, dat, unit = "trit", return = ret,
+                               weight = w, normalize = 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, normalize = TRUE)
+      set.seed(1)
+      mc <- QuartetConcordance(tree, dat, unit = "quartet", return = ret,
+                               weight = w, normalize = 8000)
+      expect_identical(is.na(exact), is.na(mc))
+      expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04)
+    }
+  }
+
+  # normalize = FALSE leaves the published raw measure untouched.
+  for (ret in c("edge", "char")) for (w in c(TRUE, FALSE)) {
+    expect_identical(
+      QuartetConcordance(tree, dat, unit = "quartet", return = ret, weight = w),
+      QuartetConcordance(tree, dat, unit = "quartet", return = ret, weight = w,
+                         normalize = FALSE))
+  }
+})
+
+test_that("QuartetConcordance() unit = 'trit' 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 = "trit",
+                              normalize = TRUE)
+  cNone <- QuartetConcordance(tree, conflict, return = "char", unit = "trit")
+  cChar <- QuartetConcordance(tree, conflict, return = "char", unit = "trit",
+                              normalize = 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)
 })