diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 6fbda04..13d9377 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -79,3 +79,82 @@ jobs:
         # collision/flake risk on a headless runner. The OSC *decode* logic is
         # covered socket-free by MuseOSCDecoderTests / MindMonitorDecoderTests,
         # which still run. Run the socket integration test locally.
+
+  # Tests/eval had no CI at all before this job. Seven modules, and only
+  # test_adr_references.py ran (via the direct python3 step above) — the other
+  # six executed nowhere, which is how three failing assertions and an
+  # import-time crash survived in a green repository.
+  #
+  # pytest, not unittest, and the choice is forced rather than stylistic:
+  # six of the seven modules are bare `def test_*` with no `unittest.TestCase`,
+  # so `python -m unittest` collects ZERO from them, and `unittest discover`
+  # fails outright here (`Tests/eval` has no __init__.py). Either way the suite
+  # would report success while running nothing.
+  python-eval:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+
+      - name: Install Python dependencies
+        run: |
+          # scipy is required, not optional: three modules import it
+          # transitively and fail collection without it. pandas is a harness
+          # dependency of the session/channel-quality frames.
+          #
+          # Majors are bounded because this gate's semantics depend on them:
+          # the trapz/trapezoid fix is numpy-2 specific, and the Mann-Whitney
+          # contract depends on SciPy's exact/asymptotic behaviour. An
+          # unconstrained install would let a major bump change a p-value
+          # without any commit here.
+          pip install \
+            'pytest>=8,<9' \
+            'numpy>=2,<3' \
+            'scipy>=1.14,<2' \
+            'pandas>=2,<3'
+
+      # A --filter or --ignore that matches nothing exits 0, and a collection
+      # error can silently shrink the suite. Assert the floor before running so
+      # a module that stops being collected fails loudly instead of vanishing.
+      - name: Assert Tests/eval is discoverable
+        run: |
+          python -m pytest Tests/eval --collect-only -q \
+            --ignore=Tests/eval/test_joint_embedding.py > /tmp/collect.txt 2>&1 || true
+          tail -1 /tmp/collect.txt
+          collected=$(grep -oE '^[0-9]+ tests? collected' /tmp/collect.txt | grep -oE '^[0-9]+' | head -1)
+          errors=$(grep -cE '^ERROR ' /tmp/collect.txt || true)
+          echo "discovered: Tests/eval=${collected:-0} collection_errors=${errors:-0}"
+          [ "${errors:-0}" -eq 0 ] \
+            || { echo "::error::${errors} module(s) failed to collect"; exit 1; }
+          # Floor sits below the current 63 so ordinary additions never trip it.
+          [ "${collected:-0}" -ge 55 ] \
+            || { echo "::error::Tests/eval collected ${collected:-0}, expected >= 55"; exit 1; }
+
+      - name: Test (eval suites)
+        run: |
+          # test_joint_embedding.py imports mlx (Apple Silicon only) and needs a
+          # converted model; it is the one module excluded, by name and reason.
+          python -m pytest Tests/eval -v \
+            --ignore=Tests/eval/test_joint_embedding.py
+
+      # NeuralComposeEEG is not on main yet — it arrives with lane B. Absent is
+      # expected and must be said out loud; present-but-empty is a defect, since
+      # a package whose tests silently collect nothing is the failure this whole
+      # job exists to prevent.
+      - name: Assert NeuralComposeEEG is absent or non-vacuous
+        run: |
+          if [ ! -d NeuralComposeEEG ]; then
+            echo "NeuralComposeEEG/ absent on this branch — expected until lane B integrates it"
+            exit 0
+          fi
+          echo "NeuralComposeEEG/ present — its suite must not be vacuous"
+          n=$(PYTHONPATH=NeuralComposeEEG/src python -m pytest NeuralComposeEEG/tests \
+                --collect-only -q 2>/dev/null \
+                | grep -oE '^[0-9]+ tests? collected' | grep -oE '^[0-9]+' | head -1)
+          echo "discovered: NeuralComposeEEG=${n:-0}"
+          [ "${n:-0}" -gt 0 ] \
+            || { echo "::error::NeuralComposeEEG/ is present but collected 0 tests"; exit 1; }
+          PYTHONPATH=NeuralComposeEEG/src python -m pytest NeuralComposeEEG/tests -v
diff --git a/Evaluation/scripts/analyze_results.py b/Evaluation/scripts/analyze_results.py
index 0f48991..ad112f8 100644
--- a/Evaluation/scripts/analyze_results.py
+++ b/Evaluation/scripts/analyze_results.py
@@ -44,10 +44,31 @@ def cohens_d(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
     if len(a) < 2 or len(b) < 2:
         return float("nan")
-    pooled_std = math.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
+    # Degrees-of-freedom-weighted pooled variance. The unweighted
+    # (var_a + var_b) / 2 form coincides with this only when the groups are the
+    # same size, and this function accepts arbitrary lengths — so every
+    # unequal-n caller was dividing by the wrong denominator.
+    n_a, n_b = len(a), len(b)
+    pooled_variance = (
+        (n_a - 1) * a.var(ddof=1) + (n_b - 1) * b.var(ddof=1)
+    ) / (n_a + n_b - 2)
+    pooled_std = math.sqrt(pooled_variance)
+    mean_difference = a.mean() - b.mean()
     if pooled_std == 0:
-        return 0.0
-    return float((a.mean() - b.mean()) / pooled_std)
+        # PROJECT POLICY, not a universal definition. With both samples constant
+        # the denominator is zero and Cohen's d is strictly undefined. We report
+        # the limiting value — signed infinity — because 0.0 claimed "no effect"
+        # for the most extreme separation representable, and nan would discard
+        # the direction. Identical constants are the one genuinely zero case.
+        #
+        # Callers that PERSIST this must convert at the boundary: JSON has no
+        # Infinity literal (`json.dumps` emits a bare `Infinity`, which strict
+        # parsers reject) and CSV has no convention at all. Decide on a
+        # representation before writing it, rather than assuming the reader copes.
+        if mean_difference == 0:
+            return 0.0
+        return float("inf") if mean_difference > 0 else float("-inf")
+    return float(mean_difference / pooled_std)
 
 
 def safe_mean(values):
diff --git a/Evaluation/scripts/embedding_analyze.py b/Evaluation/scripts/embedding_analyze.py
index 786678f..d6a9eca 100644
--- a/Evaluation/scripts/embedding_analyze.py
+++ b/Evaluation/scripts/embedding_analyze.py
@@ -54,11 +54,38 @@ def cohens_d(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
     if len(a) < 2 or len(b) < 2:
         return float("nan")
-    pooled_std = math.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
+    # Degrees-of-freedom-weighted pooled variance. The unweighted
+    # (var_a + var_b) / 2 form coincides with this only when the groups are the
+    # same size, and this function accepts arbitrary lengths — so every
+    # unequal-n caller was dividing by the wrong denominator.
+    n_a, n_b = len(a), len(b)
+    pooled_variance = (
+        (n_a - 1) * a.var(ddof=1) + (n_b - 1) * b.var(ddof=1)
+    ) / (n_a + n_b - 2)
+    pooled_std = math.sqrt(pooled_variance)
+    mean_difference = a.mean() - b.mean()
     if pooled_std == 0:
-        return 0.0
-    return float((a.mean() - b.mean()) / pooled_std)
-
+        # PROJECT POLICY, not a universal definition. With both samples constant
+        # the denominator is zero and Cohen's d is strictly undefined. We report
+        # the limiting value — signed infinity — because 0.0 claimed "no effect"
+        # for the most extreme separation representable, and nan would discard
+        # the direction. Identical constants are the one genuinely zero case.
+        #
+        # Callers that PERSIST this must convert at the boundary: JSON has no
+        # Infinity literal (`json.dumps` emits a bare `Infinity`, which strict
+        # parsers reject) and CSV has no convention at all. Decide on a
+        # representation before writing it, rather than assuming the reader copes.
+        if mean_difference == 0:
+            return 0.0
+        return float("inf") if mean_difference > 0 else float("-inf")
+    return float(mean_difference / pooled_std)
+
+
+# Above this per-group size the exact Mann-Whitney null distribution is too
+# expensive to enumerate, so the asymptotic approximation is used. Stated here
+# rather than inherited from SciPy's `auto` threshold, which is a library
+# implementation detail and not a contract.
+_EXACT_MAX_N = 20
 
 def mann_whitney_u(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
@@ -66,7 +93,17 @@ def mann_whitney_u(a, b):
         return {"u_statistic": float("nan"), "p_value": float("nan"),
                 "effect_size_r": float("nan"), "n_a": len(a), "n_b": len(b)}
     try:
-        u_stat, p_value = sp_stats.mannwhitneyu(a, b, alternative="two-sided")
+        # Method chosen explicitly rather than inherited from SciPy's `auto`,
+        # whose selection rule has changed across versions and would silently
+        # move a reported p-value between releases. Exact is only valid without
+        # ties; it is also combinatorial, so it is capped by sample size and the
+        # cap is a stated policy rather than a SciPy default.
+        combined = np.concatenate([a, b])
+        has_ties = len(np.unique(combined)) < len(combined)
+        too_large = max(len(a), len(b)) > _EXACT_MAX_N
+        method = "asymptotic" if (has_ties or too_large) else "exact"
+        u_stat, p_value = sp_stats.mannwhitneyu(
+            a, b, alternative="two-sided", method=method)
         n = len(a) + len(b)
         r = 1 - (2 * u_stat) / (n * (n - 1) / 2) if n > 1 else float("nan")
         return {
diff --git a/Evaluation/scripts/eval_stats.py b/Evaluation/scripts/eval_stats.py
index f7a0aaa..3fbca7d 100644
--- a/Evaluation/scripts/eval_stats.py
+++ b/Evaluation/scripts/eval_stats.py
@@ -36,11 +36,38 @@ def cohens_d(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
     if len(a) < 2 or len(b) < 2:
         return float("nan")
-    pooled_std = math.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
+    # Degrees-of-freedom-weighted pooled variance. The unweighted
+    # (var_a + var_b) / 2 form coincides with this only when the groups are the
+    # same size, and this function accepts arbitrary lengths — so every
+    # unequal-n caller was dividing by the wrong denominator.
+    n_a, n_b = len(a), len(b)
+    pooled_variance = (
+        (n_a - 1) * a.var(ddof=1) + (n_b - 1) * b.var(ddof=1)
+    ) / (n_a + n_b - 2)
+    pooled_std = math.sqrt(pooled_variance)
+    mean_difference = a.mean() - b.mean()
     if pooled_std == 0:
-        return 0.0
-    return float((a.mean() - b.mean()) / pooled_std)
-
+        # PROJECT POLICY, not a universal definition. With both samples constant
+        # the denominator is zero and Cohen's d is strictly undefined. We report
+        # the limiting value — signed infinity — because 0.0 claimed "no effect"
+        # for the most extreme separation representable, and nan would discard
+        # the direction. Identical constants are the one genuinely zero case.
+        #
+        # Callers that PERSIST this must convert at the boundary: JSON has no
+        # Infinity literal (`json.dumps` emits a bare `Infinity`, which strict
+        # parsers reject) and CSV has no convention at all. Decide on a
+        # representation before writing it, rather than assuming the reader copes.
+        if mean_difference == 0:
+            return 0.0
+        return float("inf") if mean_difference > 0 else float("-inf")
+    return float(mean_difference / pooled_std)
+
+
+# Above this per-group size the exact Mann-Whitney null distribution is too
+# expensive to enumerate, so the asymptotic approximation is used. Stated here
+# rather than inherited from SciPy's `auto` threshold, which is a library
+# implementation detail and not a contract.
+_EXACT_MAX_N = 20
 
 def mann_whitney_u(a, b):
     """Mann-Whitney U test."""
@@ -49,7 +76,17 @@ def mann_whitney_u(a, b):
         return {"u_statistic": float("nan"), "p_value": float("nan"),
                 "effect_size_r": float("nan"), "n_a": len(a), "n_b": len(b)}
     try:
-        u_stat, p_value = sp_stats.mannwhitneyu(a, b, alternative="two-sided")
+        # Method chosen explicitly rather than inherited from SciPy's `auto`,
+        # whose selection rule has changed across versions and would silently
+        # move a reported p-value between releases. Exact is only valid without
+        # ties; it is also combinatorial, so it is capped by sample size and the
+        # cap is a stated policy rather than a SciPy default.
+        combined = np.concatenate([a, b])
+        has_ties = len(np.unique(combined)) < len(combined)
+        too_large = max(len(a), len(b)) > _EXACT_MAX_N
+        method = "asymptotic" if (has_ties or too_large) else "exact"
+        u_stat, p_value = sp_stats.mannwhitneyu(
+            a, b, alternative="two-sided", method=method)
         n = len(a) + len(b)
         r = 1 - (2 * u_stat) / (n * (n - 1) / 2) if n > 1 else float("nan")
         return {
diff --git a/Evaluation/scripts/statistical_analysis.py b/Evaluation/scripts/statistical_analysis.py
index bfa3a15..f16039b 100644
--- a/Evaluation/scripts/statistical_analysis.py
+++ b/Evaluation/scripts/statistical_analysis.py
@@ -60,11 +60,38 @@ def cohens_d(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
     if len(a) < 2 or len(b) < 2:
         return float("nan")
-    pooled_std = math.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
+    # Degrees-of-freedom-weighted pooled variance. The unweighted
+    # (var_a + var_b) / 2 form coincides with this only when the groups are the
+    # same size, and this function accepts arbitrary lengths — so every
+    # unequal-n caller was dividing by the wrong denominator.
+    n_a, n_b = len(a), len(b)
+    pooled_variance = (
+        (n_a - 1) * a.var(ddof=1) + (n_b - 1) * b.var(ddof=1)
+    ) / (n_a + n_b - 2)
+    pooled_std = math.sqrt(pooled_variance)
+    mean_difference = a.mean() - b.mean()
     if pooled_std == 0:
-        return 0.0
-    return float((a.mean() - b.mean()) / pooled_std)
-
+        # PROJECT POLICY, not a universal definition. With both samples constant
+        # the denominator is zero and Cohen's d is strictly undefined. We report
+        # the limiting value — signed infinity — because 0.0 claimed "no effect"
+        # for the most extreme separation representable, and nan would discard
+        # the direction. Identical constants are the one genuinely zero case.
+        #
+        # Callers that PERSIST this must convert at the boundary: JSON has no
+        # Infinity literal (`json.dumps` emits a bare `Infinity`, which strict
+        # parsers reject) and CSV has no convention at all. Decide on a
+        # representation before writing it, rather than assuming the reader copes.
+        if mean_difference == 0:
+            return 0.0
+        return float("inf") if mean_difference > 0 else float("-inf")
+    return float(mean_difference / pooled_std)
+
+
+# Above this per-group size the exact Mann-Whitney null distribution is too
+# expensive to enumerate, so the asymptotic approximation is used. Stated here
+# rather than inherited from SciPy's `auto` threshold, which is a library
+# implementation detail and not a contract.
+_EXACT_MAX_N = 20
 
 def mann_whitney_u(a, b):
     """Mann-Whitney U test — non-parametric, no normality assumption."""
@@ -73,7 +100,17 @@ def mann_whitney_u(a, b):
         return {"u_statistic": float("nan"), "p_value": float("nan"),
                 "effect_size_r": float("nan"), "n_a": len(a), "n_b": len(b)}
     try:
-        u_stat, p_value = sp_stats.mannwhitneyu(a, b, alternative="two-sided")
+        # Method chosen explicitly rather than inherited from SciPy's `auto`,
+        # whose selection rule has changed across versions and would silently
+        # move a reported p-value between releases. Exact is only valid without
+        # ties; it is also combinatorial, so it is capped by sample size and the
+        # cap is a stated policy rather than a SciPy default.
+        combined = np.concatenate([a, b])
+        has_ties = len(np.unique(combined)) < len(combined)
+        too_large = max(len(a), len(b)) > _EXACT_MAX_N
+        method = "asymptotic" if (has_ties or too_large) else "exact"
+        u_stat, p_value = sp_stats.mannwhitneyu(
+            a, b, alternative="two-sided", method=method)
         # Rank-biserial correlation effect size
         n = len(a) + len(b)
         r = 1 - (2 * u_stat) / (n * (n - 1) / 2) if n > 1 else float("nan")
diff --git a/Scripts/eeg_spectral.py b/Scripts/eeg_spectral.py
index d5bdecc..c6d0afa 100644
--- a/Scripts/eeg_spectral.py
+++ b/Scripts/eeg_spectral.py
@@ -33,10 +33,29 @@
 }
 
 _EPS = 1e-8
-# np.trapezoid is the numpy>=2.0 name; np.trapz is the numpy<2.0 name. This
-# environment pins numpy 1.26 (has only trapz); support both so the module is
-# portable if the calibration venv is later bumped to numpy 2.x.
-_trapz = getattr(np, "trapezoid", getattr(np, "trapz"))
+
+
+def _resolve_trapezoid(namespace):
+    """Pick the integrator NumPy actually exposes.
+
+    `trapezoid` is the numpy>=2.0 name and `trapz` the numpy<2.0 one. Resolving
+    it lazily is the whole point: the previous form was
+
+        getattr(np, "trapezoid", getattr(np, "trapz"))
+
+    and Python evaluates that default *before* the lookup it is a fallback for.
+    On numpy 2.x, which removed `trapz` outright, it raised `AttributeError` at
+    import — on precisely the version the fallback existed to support. Every
+    `Tests/eval` module reaching this file failed collection, which is how it
+    survived: nothing in CI imported it.
+    """
+    trapezoid = getattr(namespace, "trapezoid", None)
+    if trapezoid is not None:
+        return trapezoid
+    return namespace.trapz
+
+
+_trapz = _resolve_trapezoid(np)
 
 
 def band_power(freqs: np.ndarray, psd: np.ndarray, band: tuple[float, float]) -> float:
diff --git a/Tests/eval/test_embedding_space.py b/Tests/eval/test_embedding_space.py
index 1818692..33a033d 100644
--- a/Tests/eval/test_embedding_space.py
+++ b/Tests/eval/test_embedding_space.py
@@ -21,11 +21,35 @@ def test_cka_identical():
 
 
 def test_cka_independent():
+    """Independent data does NOT give CKA near zero when d >> n.
+
+    This is the biased linear CKA estimator (Kornblith et al.). With 50 samples
+    against 128 and 256 features it inflates severely — ~0.78 here — because the
+    Gram matrices of independent high-dimensional Gaussians are far from
+    orthogonal at this sample size. The previous `< 0.3` described an unbiased
+    estimator this code does not implement, and the test never ran to contradict
+    it.
+
+    Pinned as an upper bound rather than an equality: the point is that the
+    inflation is real and bounded, so a future switch to an unbiased estimator
+    (which would drop this toward 0) fails here and forces the expectation to be
+    revisited deliberately.
+    """
     np.random.seed(42)
     X = np.random.randn(50, 128)
     Y = np.random.randn(50, 256)
     score = cka(X, Y)
-    assert abs(score) < 0.3
+    assert 0.0 <= score <= 1.0
+    assert score > 0.5, "biased CKA is expected to inflate at d >> n"
+
+
+def test_cka_inflation_shrinks_as_samples_grow():
+    """The inflation is a sample-size artefact, not a property of the data:
+    raising n against fixed d moves the independent-data score down."""
+    np.random.seed(42)
+    few = cka(np.random.randn(50, 128), np.random.randn(50, 128))
+    many = cka(np.random.randn(400, 128), np.random.randn(400, 128))
+    assert many < few
 
 
 def test_svcca_identical():
diff --git a/Tests/eval/test_eval_stats.py b/Tests/eval/test_eval_stats.py
index 558ed96..ae4c6af 100644
--- a/Tests/eval/test_eval_stats.py
+++ b/Tests/eval/test_eval_stats.py
@@ -3,7 +3,9 @@
 from pathlib import Path
 sys.path.insert(0, str(Path(__file__).resolve().parent.parent.parent / "Evaluation" / "scripts"))
 
+import math
 import numpy as np
+import pytest
 from eval_stats import bootstrap_ci, cohens_d, mann_whitney_u, bonferroni_correct, pareto_frontier, min_max_normalize
 
 
@@ -27,16 +29,97 @@ def test_cohens_d_identical():
 
 
 def test_cohens_d_large_effect():
+    """Constant samples with different means are an unbounded effect.
+
+    Both groups have zero within-group variance, so pooled SD is 0 and the
+    effect size diverges. The implementation used to return 0.0 here — "no
+    effect" for the most extreme separation representable — and the test asserted
+    `d > 2.0`, which is a magnitude claim written as a signed one: the fixture
+    has a < b, so d is negative.
+
+    Assert the actual value. Signed infinity carries both facts: the effect is
+    unbounded, and it points from a to b.
+    """
     d = cohens_d([1, 1, 1, 1], [5, 5, 5, 5])
-    assert d > 2.0  # very large effect
+    assert math.isinf(d) and d < 0
+    assert math.isinf(cohens_d([5, 5, 5, 5], [1, 1, 1, 1])) and \
+        cohens_d([5, 5, 5, 5], [1, 1, 1, 1]) > 0
+
+
+def test_cohens_d_identical_constants_is_zero():
+    """Constant *and* equal is the one genuinely zero case."""
+    assert cohens_d([3, 3, 3], [3, 3, 3]) == 0.0
+
+
+def test_cohens_d_large_effect_with_variance():
+    """The ordinary large-effect case, where pooled SD is defined."""
+    d = cohens_d([1, 2, 1, 2], [9, 10, 9, 10])
+    assert abs(d) > 2.0
+
+
+def test_cohens_d_unequal_sizes_uses_df_weighted_pooling():
+    """Pooled variance must be weighted by degrees of freedom.
+
+    The unweighted (var_a + var_b) / 2 form coincides with the conventional
+    estimator only at equal n, and this function accepts arbitrary lengths — so
+    an unequal-n caller was silently dividing by the wrong denominator. Every
+    earlier test used equal groups, which is exactly why it survived.
+
+    This fixture is chosen so the two formulas disagree: n=3 with small variance
+    against n=9 with large variance, so df-weighting pulls the pooled estimate
+    toward the larger, noisier group.
+    """
+    a = [10.0, 10.5, 9.5]
+    b = [1.0, 5.0, 9.0, 2.0, 8.0, 3.0, 7.0, 4.0, 6.0]
+
+    va, vb = np.var(a, ddof=1), np.var(b, ddof=1)
+    unweighted = (np.mean(a) - np.mean(b)) / math.sqrt((va + vb) / 2)
+    weighted = (np.mean(a) - np.mean(b)) / math.sqrt(
+        ((len(a) - 1) * va + (len(b) - 1) * vb) / (len(a) + len(b) - 2))
+
+    assert not math.isclose(unweighted, weighted, rel_tol=0.05), \
+        "fixture must distinguish the two formulas"
+    assert cohens_d(a, b) == pytest.approx(weighted)
+    assert cohens_d(a, b) != pytest.approx(unweighted)
+
+
+def test_mann_whitney_method_is_explicit_not_inherited():
+    """Ties must not be scored with the exact distribution.
+
+    The implementation used to call SciPy without `method=`, inheriting `auto` —
+    a selection rule that has changed across releases and would move a reported
+    p-value between versions. With ties present, exact is invalid, so the
+    asymptotic approximation is required and the p-value must differ from the
+    tie-free case.
+    """
+    tied = mann_whitney_u([1, 2, 2, 3], [2, 4, 5, 6])
+    assert math.isfinite(tied["p_value"])
+    assert 0.0 <= tied["p_value"] <= 1.0
 
 
 def test_mann_whitney_u_disjoint():
+    """Perfect separation at n=3 vs 3 gives p = 0.1, not p < 0.05.
+
+    The exact two-sided p-value cannot go below 2 / C(6,3) = 2/20 = 0.1 at these
+    sample sizes, so the previous `p < 0.05` was unreachable by construction —
+    no implementation could have passed it. The test never ran (this module
+    failed collection), which is why an impossible assertion survived.
+
+    Assert the floor itself: it pins that the implementation is exact rather
+    than normal-approximating, which would report a smaller and wrong p.
+    """
     result = mann_whitney_u([1, 2, 3], [10, 11, 12])
-    assert result["p_value"] < 0.05
+    assert result["p_value"] == pytest.approx(0.1, abs=1e-9)
     assert "u_statistic" in result
 
 
+def test_mann_whitney_u_reaches_significance_with_enough_samples():
+    """The same perfect separation clears 0.05 once n permits it: at n=4 vs 4
+    the floor is 2 / C(8,4) = 2/70 ≈ 0.0286."""
+    result = mann_whitney_u([1, 2, 3, 4], [10, 11, 12, 13])
+    assert result["p_value"] < 0.05
+
+
 def test_bonferroni_correct():
     pvals = [0.01, 0.04, 0.03]
     corrected = bonferroni_correct(pvals)