Anchor the expected-MI recurrence at the mode; guard expected_mi() and quartet_concordance() inputs

[ms609]

Fixes #91
Fixes #104
Fixes #105
Fixes #107

Four defects in src/expected_mi.cpp and src/quartet_concordance.cpp. No
exported signature changes: Rcpp::compileAttributes() on this tree returns
src/RcppExports.cpp and R/RcppExports.R byte-identical, so src/ts_rcpp.cpp,
src/TreeSearch-init.c and the RcppExports pair are untouched and this branch
does not collide with #84, #122 or #56.

#91 — the underflow, and why anchoring at the mode fixes it

expected_mi() sums, over each block mj of the second partition, the mutual
information contributed by an overlap k with the first partition's block a,
weighted by the hypergeometric probability

P(K = k) = C(mj, k) C(N - mj, a - k) / C(N, a).

The old code evaluated P(K = kmin) once by exponentiating a sum of log
factorials, then walked upwards with the exact ratio
P(k+1)/P(k) = (mj - k)(a - k) / ((k + 1)(N - mj - a + k + 1)).

The seed is the problem. kmin = max(0, a + mj - N), so for a near-balanced
split against a near-balanced character kmin = 0 and P(K = 0) = 1 / C(N, a),
which is about 2^-1197 at N = 1200. That is below the smallest subnormal
double (2^-1074), so it becomes exactly zero; the recurrence is
multiplicative, so every subsequent term is zero too, and the if (Pk > 0.0)
guard suppresses that block's whole contribution.

The symptom is worse than "returns 0". The seed is evaluated per block mj,
so a character whose blocks straddle the threshold has some blocks silently
dropped from the sum and others counted, and expected_mi() returns a
plausible-looking but truncated number. Over 600 random partitions with N from
10 to 3000 (below), 60 came back wrong; of those only 4 were exactly zero, and
the rest ranged from a part in a million down to 26 % of the true value.
ClusteringConcordance(normalize = TRUE) subtracted whatever survived.

The fix keeps the same recurrence but seeds it at the mode
kmode = floor((mj + 1)(a + 1) / (N + 2)), clamped to [kmin, kmax], and walks
outwards in both directions (the downward step uses the reciprocal ratio,
P(k-1)/P(k) = k(N - mj - a + k) / ((mj - k + 1)(a - k + 1))). P(K = kmode)
is the largest of at most N + 1 probabilities that sum to one, so it is at
least 1/(N + 1) and can never underflow — the seed is now unconditionally
representable rather than accidentally so. The total number of k values
visited is unchanged, so the per-k cost the recurrence was introduced for
(commit 7904371, 2025-09-01) is preserved in full. Because the distribution is
log-concave and kmode is an argmax, once a tail term reaches exactly zero
every remaining term in that direction is also exactly zero under the same
multiplicative recurrence, so each sweep breaks there; the break is exactly
equivalent to running on, not merely a good approximation.

One incidental hazard goes with it: the ratio was formed as
static_cast<double>((mj - k) * (a - k)), i.e. the product was taken in int
and overflows above roughly a·mj > 2^31; it is now formed in double.

Evidence

Two reference implementations: ReferenceEmi() computes the same sum from R's
lchoose(), which works in log space and cannot underflow at the tails, and a
second check recomputes it from dhyper(). Both are independent implementations
but share the closed form, so they test the recurrence's numerics, not the
formula. The formula itself is checked separately, by the asymptote below and by
a formula-free permutation estimate.

Correctness above the old ceiling. For two balanced binary partitions the
expected mutual information has the closed-form asymptote 1/(2N ln 2) bits
(2N·MI in nats converges on a 1-df chi-square, whose mean is 1), which the
computed values approach as 1 + 1.5/N:

N	expected_mi	1/(2N ln 2)	ratio
100	0.00732365294	0.00721347520	1.01527
500	0.00144703860	0.00144269504	1.00301
1000	0.00072243147	0.00072134752	1.00150
1200	0.00060187545	0.00060112293	1.00125
2000	0.00036094451	0.00036067376	1.00075
5000	0.00014431280	0.00014426950	1.00030

This is the independent test of the form, and it is the strongest evidence
here: the fixed values are not merely non-zero, they land on an analytic
prediction derived without reference to the code. They are also strictly
decreasing in N across 100…5000, as a chance correction must be.

Agreement with the log-space oracle.

• 20 540 exhaustive 2×2 cases (every N in 2…40, every a, every mj):
  worst relative deviation 2.8e-13. A k counted twice or missed would show up
  here immediately.
• 600 random partitions, N in {10, 30, 100, 300, 700, 1000, 1200, 1500, 2000,
  3000}, nj of 2–6 blocks: worst relative deviation 3.6e-11 (the largest at
  N = 3000, where the accumulated error of the 8192-entry log-factorial table
  dominates).
• 4000 random partitions, N uniform on 4…900: worst 4.7e-12.

Agreement with the pre-fix values where the pre-fix code was right. Pre-fix
values come from origin/cpp-search's sources compiled unchanged into a
separate library. Of the 600 random partitions, 540 were accurate pre-fix
(within 1e-10 of the oracle); on those the fixed value differs from the pre-fix
value by at most 1.5e-11 relative, and the two bracket the oracle rather
than one being systematically nearer it (worst |pre-fix − oracle| 2.8e-11,
worst |fixed − oracle| 2.0e-11, neither consistently smaller). Where the old
code was right, the new code is right to the same order — this is not a change
of answers, it is the same recurrence started somewhere it can be started.
These are sample maxima, not proven bounds; an independent replication by a
reviewer at N ≤ 900 measured 2.3e-12 against my 2.2e-12 on the same regime.

Worked examples:

case	pre-fix	fixed	oracle
ni={3,4} nj={2,5} N=7	0.15383715015513183	0.15383715015513186	0.15383715015513161
ni={50,50} nj={50,50} N=100	0.007323652940324857	0.0073236529403252724	0.0073236529403244337
ni={500,500} nj={500,500} N=1000	0.00072243147032421289	0.00072243147032464647	0.00072243147032802376
ni={400,800} nj={300,400,500} N=1200	0.0012047105794341356	0.0012047105794329065	0.0012047105794319285
ni={550,550} nj={550,550} N=1100	0	0.00065666615423653921	0.00065666615423918163
ni={600,600} nj={600,600} N=1200	0	0.00060187545297477572	0.00060187545297570337
ni={1000,1000} nj={1000,1000} N=2000	0	0.00036094450630637685	0.00036094450630407357
ni={4500,4500} nj={4500,4500} N=9000	0	8.0163085423003874e-05	8.0163085420675143e-05

Who was affected. The relevant N is the number of tips a character
scores, not tips in the tree (R/Concordance.R:216, tabulate() over non-NA
states), so heavy missing data lowers the exposure. Against a perfectly even
character the first spurious zero appears at N = 1080; sweeping every split
size a:

N	split sizes wrongly zero	affected range of a
500	0 / 499	—
1000	0 / 999	—
1100	5 / 1099	548–552
1200	39 / 1199	581–619
1500	187 / 1499	657–843
2000	519 / 1999	741–1259

The band is centred on the balanced splits — the deep edges of a large tree,
where the correction matters most — and it widens with N. Against uneven
characters the truncation is partial rather than total, so the count of exact
zeros understates the reach: on the 600-partition sample, 1 of 60 partitions was
wrong at N = 1200, 13 of 60 at 1500, 30 of 60 at 3000.

#104 — ni[1] read out of bounds

expected_mi() now stops with "ni must be a vector of length 2." when
ni.size() != 2, matching the guard mi_key() already had and the contract
man/expected_mi.Rd has always documented ("Integer vector of length two").
This is not theoretical: on the pre-fix build, expected_mi(integer(0), c(2L, 5L))
reads two garbage ints and then indexes the log-factorial table with them —
at -O2 -DNDEBUG, the pre-fix test run segfaults at that line (exit 139).
#104 and #105 compose into a crash; either guard alone stops it.

#105 — portability and the unguarded table index

__attribute__((constructor)) is a GCC/Clang extension. It is replaced by a
block-scope static const std::vector<double> initialised by a lambda on first
use. DESCRIPTION declares SystemRequirements: C++17, and C++11 onwards
guarantee that such an initialisation runs exactly once even if several threads
enter concurrently, the others blocking until it completes — so this is not a
thread-safety regression relative to the constructor attribute, which ran before
any thread existed. The table is read-only afterwards, so concurrent
expected_mi() calls remain safe. (No caller runs it concurrently today; the
guarantee matters only if one later does.) l2factorial() now rejects a
negative argument rather than indexing the table below zero.

#107 — negative state codes

quartet_concordance() uses each non-NA state code directly as an index into
its n0/n1 count buffers. It now stops with "characters must contain
non-negative state codes." rather than trusting the R-side convention. The
guard sits after the existing is_na test, so NA_integer_ (INT_MIN) is
still treated as a missing entry, not a negative code. Pre-fix, a -1 in the
matrix silently returned all-zero counts (an out-of-bounds read/write,
undefined behaviour that happened not to fault here).

Tests

New file tests/testthat/test-expected-mi.R, 28 assertions, Tier 1 (arithmetic
only, whole file under a second). Verified against a build of
origin/cpp-search's unmodified sources:

test	pre-fix result
large balanced partitions correct (#91)	FAIL — 0 at N = 1100, 1200, 2000, and through .ExpectedMI(); not all positive; not decreasing
small partitions unchanged (#91)	passes pre-fix by construction — the no-regression half of #91's test, pinning the pre-fix numbers to 1e-10
ni must be a pair (#104)	FAIL — no error raised; integer(0) then segfaults the session
factorial lookup boundary, N = 9000 (#105)	FAIL — returns 0
negative state codes rejected (#107)	FAIL — no error raised, all-zero counts returned

Two caveats stated plainly rather than papered over. First, the small-partition
test cannot fail before the fix: its whole purpose is to pin the pre-fix
numbers. Second, #105's substance is a compile-time portability change plus a
guard on a branch that is unreachable through expected_mi()'s own arithmetic
once #104's length check is in place; its test exercises the replacement
initialiser across both the table and lgamma() branches of l2factorial(),
and fails pre-fix for #91's reason rather than #105's.

Local: test-expected-mi.R 28/28 pass, test-Concordance.R 75/75 pass (2 CRAN
skips), against a tarball build installed to a per-agent library.
spelling::spell_check_package(vignettes = TRUE) clean. No test uses threads or
costs more than a second.

Reviewed, not fixed here

mi_key() narrows block sizes to uint16_t, so .ExpectedMICache keys alias
above 65 535 items. Pre-existing, out of reach of the tree sizes this patch
targets, and touching it would change a signature — left alone deliberately.