diff --git a/docs/api.rst b/docs/api.rst
index 8580cb8..43ba3ab 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -115,6 +115,7 @@ API
    stats.ensure_2d
    stats.q_profile
    stats.q_gen
+   stats.log_chi2_sf
    stats.bonferroni
    stats.fdr
 
diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py
index 0a540e5..93f4af1 100644
--- a/pymare/estimators/combination.py
+++ b/pymare/estimators/combination.py
@@ -5,11 +5,10 @@
 from abc import abstractmethod
 
 import numpy as np
-import scipy.stats as ss
-from scipy.special import log_ndtr, ndtr
+from scipy.special import log_ndtr, ndtri_exp
 
 from ..results import CombinationTestResults
-from ..stats import encode_groups, normalize_group_weights
+from ..stats import encode_groups, log_chi2_sf, normalize_group_weights
 from .estimators import BaseEstimator
 
 
@@ -85,12 +84,17 @@ def __init__(self, mode="directed"):
         self.mode = mode
 
     @abstractmethod
-    def p_value(self, z, *args, **kwargs):
-        """Calculate p-values."""
+    def log_p_value(self, z, *args, **kwargs):
+        """Calculate natural logarithms of the p-values."""
         pass
 
-    def _z_to_p(self, z):
-        return ndtr(-z)
+    def p_value(self, z, *args, **kwargs):
+        """Calculate p-values.
+
+        Underflows to zero where the combined evidence exceeds what a double can
+        represent; :meth:`log_p_value` is the same quantity without that limit.
+        """
+        return np.exp(self.log_p_value(z, *args, **kwargs))
 
     def fit(self, z, *args, **kwargs):
         """Fit the estimator to z-values."""
@@ -102,18 +106,23 @@ def fit(self, z, *args, **kwargs):
             # aggregation) while evaluating the two directed tails.
             ose = copy.copy(self)
             ose.mode = "directed"
-            p1 = ose.p_value(z, *args, **kwargs)
-            p2 = ose.p_value(-z, *args, **kwargs)
-            p = np.minimum(1, 2 * np.minimum(p1, p2))
-            z_calc = ss.norm.isf(p)
-            z_calc[p2 < p1] *= -1
+            log_p1 = ose.log_p_value(z, *args, **kwargs)
+            log_p2 = ose.log_p_value(-z, *args, **kwargs)
+            # Doubling the smaller tail and capping at 1, in logs: the
+            # correction for two tests is an added log(2), the cap a minimum
+            # against log(1).
+            log_p = np.minimum(0.0, np.log(2.0) + np.minimum(log_p1, log_p2))
+            z_calc = -ndtri_exp(log_p)
+            z_calc[log_p2 < log_p1] *= -1
         else:
             if self.mode == "undirected":
                 z = np.abs(z)
-            p = self.p_value(z, *args, **kwargs)
-            z_calc = ss.norm.isf(p)
+            log_p = self.log_p_value(z, *args, **kwargs)
+            # ``norm.isf(p)`` instead would saturate to +/-inf the moment p
+            # underflowed or hit 1; the inverse of the log CDF does not.
+            z_calc = -ndtri_exp(log_p)
 
-        self.params_ = {"p": p, "z": z_calc}
+        self.params_ = {"p": np.exp(log_p), "logp": log_p, "z": z_calc}
         return self
 
     def summary(self):
@@ -124,8 +133,12 @@ def summary(self):
                 "This {} instance hasn't been fitted yet. Please "
                 "call fit() before summary().".format(name)
             )
+        # p is exactly ``exp(logp)``, which the container derives itself, so
+        # passing it would store a second copy of one quantity. z is not
+        # derivable here: in concordant mode it carries the sign of whichever
+        # tail won, which the log p-value alone does not record.
         return CombinationTestResults(
-            self, self.dataset_, z=self.params_["z"], p=self.params_["p"]
+            self, self.dataset_, z=self.params_["z"], logp=self.params_["logp"]
         )
 
 
@@ -345,8 +358,8 @@ def fit(self, z, w=None, g=None, corr=None):
         self.corr_ = corr
         return super().fit(z, w=w, g=g, corr=corr)
 
-    def p_value(self, z, w=None, g=None, corr=None):
-        """Calculate p-values."""
+    def log_p_value(self, z, w=None, g=None, corr=None):
+        """Calculate natural logarithms of the p-values."""
         if w is None:
             w = np.ones_like(z)
         else:
@@ -370,7 +383,7 @@ def p_value(self, z, w=None, g=None, corr=None):
             group_z, group_w = self._group_statistics(z, w, g, corr=corr)
             variance = np.square(group_w).sum(axis=0)
             cz = (group_z * group_w).sum(axis=0) / np.sqrt(variance)
-            return ss.norm.sf(cz)
+            return log_ndtr(-cz)
 
         if g is None and corr is not None:
             warnings.warn("Correlation matrix provided without groups. Ignoring.")
@@ -385,7 +398,10 @@ def p_value(self, z, w=None, g=None, corr=None):
         variance = (w**2).sum(0) + sigma
 
         cz = (z * w).sum(0) / np.sqrt(variance)
-        return ss.norm.sf(cz)
+        # log_ndtr, not norm.sf: the combined z is a weighted *sum*, so it grows
+        # with the number of observations and passes 38 -- where a double-
+        # precision p-value is exactly zero -- on datasets of very ordinary size.
+        return log_ndtr(-cz)
 
 
 class FisherCombinationTest(CombinationTest):
@@ -649,15 +665,13 @@ def fit(self, z, g=None, corr=None, w=None):
         self.corr_ = corr
         return super().fit(z, g=g, corr=corr, w=w)
 
-    def p_value(self, z, g=None, corr=None, w=None):
-        """Calculate p-values."""
+    def log_p_value(self, z, g=None, corr=None, w=None):
+        """Calculate natural logarithms of the p-values."""
         g, corr = self._validate_dependence_inputs(z, g, corr)
 
-        # Work in log space throughout. Going via p underflows to exactly 0
-        # around z = 38, after which log(p) is -inf and the combined result
-        # collapses to p = 0 with z = inf, no matter how many other inputs
-        # argue otherwise. log_ndtr is accurate far into that tail, so a single
-        # extreme z no longer destroys the statistic.
+        # Work in log space throughout. Going
+        # via p underflows to exactly 0 around z = 38, after which log(p) is
+        # -inf and the combined result collapses to p = 0 with z = inf.
         log_p = log_ndtr(-z)
         weights = self._group_weights(g, z.shape[0], w=w)
         if g is None and w is None:
@@ -682,4 +696,4 @@ def p_value(self, z, g=None, corr=None, w=None):
         scale = variance / (2.0 * expectation)
         dof = 2.0 * expectation**2 / variance
 
-        return ss.chi2.sf(chi2 / scale, dof)
+        return log_chi2_sf(chi2 / scale, dof)
diff --git a/pymare/results.py b/pymare/results.py
index d4d88d6..b96767a 100644
--- a/pymare/results.py
+++ b/pymare/results.py
@@ -10,6 +10,7 @@
 import numpy as np
 import pandas as pd
 import scipy.stats as ss
+from scipy.special import log_ndtr, ndtri_exp
 
 try:
     import arviz as az
@@ -25,6 +26,7 @@
     collapse_groups,
     collapse_groups_by_n,
     encode_groups,
+    log_chi2_sf,
     q_gen,
     q_profile,
 )
@@ -273,7 +275,8 @@ def get_fe_stats(self, alpha=0.05):
             est         The parameter estimate for the regressor.
             se          The standard error of the estimate.
             z           The z score of the estimate.
-            p           The p value the estimate.
+            p           The two-tailed p value of the estimate.
+            logp        The natural logarithm of ``p``.
             ci_l/ci_u   Lower and upper bounds of the estimate.
             =========== ==========================================================================
 
@@ -296,41 +299,52 @@ def get_fe_stats(self, alpha=0.05):
         how uncertain the coefficients are. The ``undefined`` branch below turns
         that into a NaN p-value rather than a maximally significant one.
 
+        The tail is computed in logs because ``p`` cannot express the
+        transform's own input: a two-tailed normal p is exactly zero from
+        ``|z| = 38.5`` on, and every deviate past that would come back
+        infinite. Via ``logp``, ``z`` grows like ``sqrt(2 |logp|)`` without
+        bound.
+
+        .. versionchanged:: 0.0.11
+            ``p`` is no longer floored at ``numpy.finfo(float).eps``. The floor
+            read a limit of precision as one of magnitude, truncating ``z`` at
+            ``norm.isf(eps / 2) = 8.21`` for every estimate below 2.2e-16,
+            however far below it fell.
+
         References
         ----------
         .. footbibliography::
 
         """
         beta, se = self.fe_params, self.fe_se
-        epsilon = np.finfo(beta.dtype).eps
-        z = beta / se
 
         # Cluster-robust standard errors are asymptotic in the number of
         # groups, so refer them to a t distribution rather than a normal.
         dof = self.fe_dof
-        if dof is None:
-            crit = ss.norm.ppf(1 - alpha / 2)
-            p = 1 - np.abs(0.5 - ss.norm.cdf(z)) * 2
-        else:
-            crit = ss.t.ppf(1 - alpha / 2, dof)
-            p = 2 * ss.t.sf(np.abs(z), dof)
+        # A zero standard error divides to +-inf, or to NaN if the estimate is
+        # zero too. Both are answered by the ``undefined`` mask below rather
+        # than by a warning, so the arithmetic is allowed to produce them.
+        with np.errstate(invalid="ignore", divide="ignore"):
+            stat = beta / se
+            if dof is None:
+                crit = ss.norm.ppf(1 - alpha / 2)
+                log_tail = ss.norm.logsf(np.abs(stat))
+                z = stat
+            else:
+                crit = ss.t.ppf(1 - alpha / 2, dof)
+                log_tail = ss.t.logsf(np.abs(stat), dof)
+                z = -np.sign(stat) * ndtri_exp(log_tail)
 
-        p = np.asarray(p, dtype=float)
-        p[p == 0] += epsilon
+            # Two tails, so an added log(2). The cap is against log(1), which
+            # the sum can exceed only by rounding at a statistic of essentially
+            # zero; it is not a floor on the tail.
+            logp = np.minimum(log_tail + math.log(2.0), 0.0)
+            p = np.exp(logp)
 
         # A standard error that is zero or non-finite carries no information,
         # but "est / 0" is +-inf and yields p = 0 -- maximally significant. Mark
         # those entries undefined instead of maximally certain.
         undefined = ~np.isfinite(se) | (se <= 0)
-        if np.any(undefined):
-            p = np.where(undefined, np.nan, p)
-
-        if dof is not None:
-            # ``est / se`` is a t statistic here, so reporting it as "z" would
-            # leave the two entries disagreeing: thresholding on the z values
-            # and on the p values would select different results. Report the
-            # z that carries the same tail probability instead.
-            z = np.sign(z) * ss.norm.isf(np.clip(p, epsilon, 1.0) / 2)
 
         stats = {
             "est": beta,
@@ -338,7 +352,8 @@ def get_fe_stats(self, alpha=0.05):
             "ci_l": np.where(undefined, np.nan, beta - crit * se),
             "ci_u": np.where(undefined, np.nan, beta + crit * se),
             "z": np.where(undefined, np.nan, z),
-            "p": p,
+            "p": np.where(undefined, np.nan, p),
+            "logp": np.where(undefined, np.nan, logp),
         }
 
         return stats
@@ -491,6 +506,8 @@ def get_heterogeneity_stats(self):
                     freedom, where n is the number of independent observations and k is the
                     number of regressors.
             p(Q)    P values associated with the Cochran's Q values.
+            logp(Q) Natural logarithms of the ``p(Q)`` values, computed directly and
+                    therefore still finite where ``p(Q)`` has underflowed to zero.
             I^2     The proportion of the variance in input estimates that is due to heterogeneity
                     instead of sampling error :footcite:p:`higgins2002quantifying`.
                     This measure is bounded from 0 to 100.
@@ -518,12 +535,18 @@ def get_heterogeneity_stats(self):
             # zero and I^2 / H are 100% / inf from pure rounding noise, which
             # reads as total heterogeneity rather than no information.
             nan = np.full(np.shape(q_fe), np.nan)
-            return {"Q": nan, "p(Q)": nan.copy(), "I^2": nan.copy(), "H": nan.copy()}
+            return {
+                "Q": nan,
+                "p(Q)": nan.copy(),
+                "logp(Q)": nan.copy(),
+                "I^2": nan.copy(),
+                "H": nan.copy(),
+            }
 
         i2 = np.maximum(100.0 * (q_fe - df) / q_fe, 0.0)
         h = np.maximum(np.sqrt(q_fe / df), 1.0)
-        p = ss.chi2.sf(q_fe, df)
-        return {"Q": q_fe, "p(Q)": p, "I^2": i2, "H": h}
+        logp = log_chi2_sf(q_fe, df)
+        return {"Q": q_fe, "p(Q)": np.exp(logp), "logp(Q)": logp, "I^2": i2, "H": h}
 
     def to_df(self, alpha=0.05):
         """Return a pandas DataFrame summarizing fixed effect results.
@@ -554,6 +577,7 @@ def to_df(self, alpha=0.05):
             se          The standard error of the estimate.
             z-score     The z score of the estimate.
             p-value     The p value the estimate.
+            -log10(p)   The p value on a base-10 log scale.
             ci_+        Lower and upper bounds of the estimate. There will be two columns, with
                         names based on the ``alpha`` value. For example, if ``alpha = 0.05``,
                         the CI columns will be ``"ci_0.025"`` and ``"ci_0.975"``.
@@ -573,10 +597,11 @@ def to_df(self, alpha=0.05):
         fe_stats = self.get_fe_stats(alpha).items()
         df = pd.DataFrame({k: v.ravel() for k, v in fe_stats})
         df["name"] = self.dataset.X_names
-        df = df.loc[:, ["name", "est", "se", "z", "p", "ci_l", "ci_u"]]
+        df["logp"] = np.abs(df["logp"]) / math.log(10.0)
+        df = df.loc[:, ["name", "est", "se", "z", "p", "logp", "ci_l", "ci_u"]]
         ci_l = "ci_{:.6g}".format(alpha / 2)
         ci_u = "ci_{:.6g}".format(1 - alpha / 2)
-        df.columns = ["name", "estimate", "se", "z-score", "p-value", ci_l, ci_u]
+        df.columns = ["name", "estimate", "se", "z-score", "p-value", "-log10(p)", ci_l, ci_u]
         return df
 
     def permutation_test(self, n_perm=1000):
@@ -752,22 +777,34 @@ class CombinationTestResults:
         Array of z-scores. Default = None.
     p : :obj:`numpy.ndarray`, optional
         Array of right-tailed p-values. Default = None.
+    logp : :obj:`numpy.ndarray`, optional
+        Array of natural logarithms of the right-tailed p-values. Default = None.
     """
 
-    def __init__(self, estimator, dataset, z=None, p=None):
+    def __init__(self, estimator, dataset, z=None, p=None, logp=None):
         self.estimator = estimator
         self.dataset = dataset
-        if p is None and z is None:
-            raise ValueError("One of 'z' or 'p' must be provided.")
+        if p is None and z is None and logp is None:
+            raise ValueError("One of 'z', 'p' or 'logp' must be provided.")
         self._z = z
         self._p = p
+        self._logp = logp
+
+    @property
+    @lru_cache(maxsize=1)
+    def logp(self):
+        """Natural logarithms of the p-values."""
+        if self._logp is None:
+            with np.errstate(divide="ignore"):
+                self._logp = log_ndtr(-self._z) if self._z is not None else np.log(self._p)
+        return self._logp
 
     @property
     @lru_cache(maxsize=1)
     def z(self):
         """Z-values."""
         if self._z is None:
-            self._z = ss.norm.isf(self.p)
+            self._z = -ndtri_exp(self.logp)
         return self._z
 
     @property
@@ -775,7 +812,7 @@ def z(self):
     def p(self):
         """P-values."""
         if self._p is None:
-            self._p = ss.norm.sf(self.z)
+            self._p = np.exp(self.logp)
         return self._p
 
     def permutation_test(self, n_perm=1000):
@@ -904,14 +941,8 @@ def permutation_test(self, n_perm=1000):
                 kwargs["corr"] = permutation_corr
             params = est.fit(**kwargs).params_
 
-            # Compare on p-values, not on z. The reported z is
-            # ``norm.isf(p)``, which saturates to +/-inf as soon as p hits 0
-            # or 1 -- routine for concordant mode, where p is capped at 1.
-            # Every permutation then ties at -inf and the comparison degrades
-            # to "always significant". p is the actual test statistic in every
-            # mode, is monotone in the evidence, and never overflows.
-            observed = np.ravel(self.p)[i]
-            null = np.ravel(params["p"])
+            observed = np.ravel(self.logp)[i]
+            null = np.ravel(params["logp"])
             p_p[i] = (null <= observed).mean()
 
         # p-values can't be smaller than 1/n_perm
diff --git a/pymare/stats.py b/pymare/stats.py
index 173f720..40e769e 100644
--- a/pymare/stats.py
+++ b/pymare/stats.py
@@ -20,6 +20,7 @@
 import numpy as np
 import scipy.stats as ss
 from scipy.optimize import Bounds, minimize
+from scipy.special import gammaln
 
 # At or below this many clusters, robust variance estimation is known to be
 # anti-conservative; see Hedges, Tipton & Johnson (2010) and Tipton (2015).
@@ -2241,6 +2242,107 @@ def q_profile(y, v, X, alpha=0.05, groups=None):
     return {"ci_l": lb, "ci_u": ub}
 
 
+#: Iterations allowed in the continued fraction of :func:`log_chi2_sf`. A safety
+#: stop rather than a working limit. Convergence slows as ``x`` approaches ``a``
+#: from above -- 722 iterations at ``a = 5e5``, ``x = a + 1`` -- but that is the
+#: regime where the tail is around one half and SciPy is exact, so the fraction
+#: is never used there. Where it *is* used, past the point at which a
+#: double-precision tail underflows, four to six iterations suffice.
+_GAMMA_CF_MAX_ITER = 300
+
+#: Relative change below which the continued fraction is treated as converged.
+#: Set at machine epsilon because each iteration multiplies the running value by
+#: a factor approaching one; stopping earlier costs digits in the log, which is
+#: the quantity this function exists to get right.
+_GAMMA_CF_TOL = np.finfo(np.float64).eps
+
+
+def log_chi2_sf(q, df):
+    r"""Natural logarithm of the chi-squared upper tail, without underflowing.
+
+    Parameters
+    ----------
+    q : :obj:`numpy.ndarray` or :obj:`float`
+        Statistic values, non-negative.
+    df : :obj:`numpy.ndarray` or :obj:`float`
+        Degrees of freedom, positive. Broadcast against ``q``.
+
+    Returns
+    -------
+    :obj:`numpy.ndarray`
+        ``log(P(chi^2_df > q))``.
+
+    Notes
+    -----
+    ``scipy.stats.chi2.logsf`` computes the survival function first and takes
+    its logarithm afterwards, so it returns ``-inf`` for every ``q`` whose tail
+    falls below the smallest positive double -- which Cochran's Q reaches on a
+    few hundred heterogeneous observations, well inside the range of an ordinary
+    meta-analysis.
+
+    The tail is the regularized upper incomplete gamma function
+    :math:`Q(a, x)` at :math:`a = df/2`, :math:`x = q/2`, which factors as
+
+    .. math::
+
+        Q(a, x) = \frac{e^{-x} x^{a}}{\Gamma(a)} \cdot F(a, x),
+
+    where :math:`F` is a continued fraction of order one, evaluated here by the
+    modified Lentz algorithm. Everything that underflows lives in the prefactor,
+    and the prefactor is an exponential, so evaluating the fraction on its own
+    and adding ``-x + a log x - lnGamma(a)`` to its logarithm gives the answer
+    with no intermediate that can flush to zero.
+
+    The fraction is used only where SciPy has already failed, and SciPy is used
+    everywhere else. That is not a preference for one over the other but the
+    split that makes both fast: the fraction's convergence slows as ``x``
+    approaches ``a`` from above, needing hundreds of iterations at ``x = a + 1``
+    with a large ``a``, and that is precisely the region where the tail is
+    around one half, nothing underflows and SciPy is exact. Where the tail is
+    small enough to have underflowed, ``x`` is far enough beyond ``a`` that the
+    fraction converges in a handful of iterations.
+    """
+    q = np.asarray(q, dtype=float)
+    df = np.asarray(df, dtype=float)
+    a, x = np.broadcast_arrays(df / 2.0, q / 2.0)
+
+    with np.errstate(divide="ignore"):
+        out = np.broadcast_to(ss.chi2.logsf(q, df), a.shape).astype(float, copy=True)
+
+    # ``x >= a + 1`` is the fraction's domain and is implied by an underflowed
+    # tail, but it is asserted rather than assumed: a NaN or infinite input
+    # reaches SciPy's answer through the same non-finite test.
+    tail = ~np.isfinite(out) & np.isfinite(x) & (x >= a + 1.0)
+    if not np.any(tail):
+        return out
+
+    a_t, x_t = a[tail], x[tail]
+    tiny = np.finfo(np.float64).tiny
+
+    # Modified Lentz: b, c and d track the fraction's recurrence, h its value.
+    b = x_t + 1.0 - a_t
+    c = np.full(b.shape, 1.0 / tiny)
+    d = 1.0 / b
+    h = d.copy()
+    active = np.ones(b.shape, dtype=bool)
+    for i in range(1, _GAMMA_CF_MAX_ITER + 1):
+        an = -i * (i - a_t)
+        b = b + 2.0
+        d = an * d + b
+        d = np.where(np.abs(d) < tiny, tiny, d)
+        c = b + an / c
+        c = np.where(np.abs(c) < tiny, tiny, c)
+        d = 1.0 / d
+        delta = d * c
+        h = np.where(active, h * delta, h)
+        active &= np.abs(delta - 1.0) > _GAMMA_CF_TOL
+        if not np.any(active):
+            break
+
+    out[tail] = -x_t + a_t * np.log(x_t) - gammaln(a_t) + np.log(h)
+    return out
+
+
 def q_gen(y, v, X, tau2, groups=None):
     """Calculate a generalized form of Cochran's Q-statistic.
 
diff --git a/pymare/tests/conftest.py b/pymare/tests/conftest.py
index 440f55e..2803055 100644
--- a/pymare/tests/conftest.py
+++ b/pymare/tests/conftest.py
@@ -60,6 +60,26 @@ def variables():
     return (y, v, X)
 
 
+@pytest.fixture(scope="package")
+def extreme_effect_results():
+    """Return the results of a fit at a chosen effect size.
+
+    The sampling variances are small and the residuals around the effect are
+    tiny, so the statistic grows in proportion to the effect and can be driven
+    as far into the tail as a caller needs -- including past the point where a
+    double-precision p-value is exactly zero.
+    """
+    v = np.array([[0.02, 0.03, 0.025, 0.04, 0.02]]).T
+    residual = np.array([[0.01, -0.02, 0.015, -0.01, 0.005]]).T
+
+    def fit(effect, correction="knapp-hartung"):
+        dataset = Dataset(y=effect + residual, v=v)
+        estimator = WeightedLeastSquares(tau2=0.0, small_sample_correction=correction)
+        return estimator.fit_dataset(dataset).summary()
+
+    return fit
+
+
 @pytest.fixture(scope="package")
 def small_variance_variables(variables):
     """Make highly correlated variables."""
diff --git a/pymare/tests/test_combination_tests.py b/pymare/tests/test_combination_tests.py
index 4e0a133..883e17e 100644
--- a/pymare/tests/test_combination_tests.py
+++ b/pymare/tests/test_combination_tests.py
@@ -348,6 +348,49 @@ def test_fisher_does_not_underflow_on_extreme_z():
     assert np.allclose(fitted["p"], 1.035e-244, rtol=1e-3)
 
 
+def test_fisher_does_not_underflow_on_many_moderate_z():
+    """Its own combined tail underflows too, and much sooner than one input can.
+
+    Converting each input in logs, as the test above checks, still leaves
+    ``chi2.sf`` -- and ``chi2.logsf``, which logs it afterwards -- at the end of
+    the statistic. Two hundred inputs at z = 3 put the combined chi-squared
+    where a double-precision tail is exactly zero, which came back as
+    ``logp = -inf`` and ``z = inf``: less informative than any single input.
+    """
+    fitted = FisherCombinationTest().fit(np.full((200, 2), 3.0)).params_
+
+    assert np.all(fitted["p"] == 0.0)  # the representation, not the evidence
+    assert np.allclose(fitted["logp"], -749.1910315, rtol=1e-8)
+    assert np.allclose(fitted["z"], 38.5906313, rtol=1e-7)
+
+
+def test_stouffer_does_not_underflow_on_many_moderate_z():
+    """Fisher had this guard already; Stouffer reaches the same wall sooner.
+
+    The combined statistic is a weighted *sum* divided by the square root of the
+    weight, so it grows like ``sqrt(k)``: four hundred inputs at z = 3 -- not an
+    unusual meta-analysis -- combine to z = 60, whose one-tailed p-value is
+    about 1e-785. Going via ``norm.sf`` returned exactly 0 there, and the z
+    rebuilt from that 0 was ``+inf``, so the answer was less informative than
+    any single one of its inputs.
+    """
+    z = np.full((400, 2), 3.0)
+
+    fitted = StoufferCombinationTest().fit(z).params_
+
+    assert np.all(fitted["p"] == 0.0)  # the representation, not the evidence
+    assert np.all(np.isfinite(fitted["logp"]))
+    assert np.allclose(fitted["z"], 60.0)
+    assert np.allclose(fitted["logp"], ss.norm.logsf(60.0))
+
+
+def test_public_p_value_survives_the_move_to_log_space():
+    """Subclasses now implement log_p_value, but p_value stays part of the API."""
+    est = StoufferCombinationTest()
+
+    assert np.allclose(est.p_value(_z2), np.exp(est.log_p_value(_z2)))
+
+
 def test_constant_estimates_cannot_yield_a_correlation(combination_estimator):
     """A row that never varies has no correlation, so do not return NaN."""
     z = np.full((4, 5), 1.3)
diff --git a/pymare/tests/test_core.py b/pymare/tests/test_core.py
index d8b99c4..099b45a 100644
--- a/pymare/tests/test_core.py
+++ b/pymare/tests/test_core.py
@@ -112,7 +112,7 @@ def test_meta_regression_2(dataset_n):
     """Test meta_regression function."""
     y, n = dataset_n.y, dataset_n.n
     df = meta_regression(y=y, n=n).to_df()
-    assert df.shape == (1, 7)
+    assert df.shape == (1, 8)
 
 
 # -----------------------------------------------------------------------------
diff --git a/pymare/tests/test_metafor_alignment.py b/pymare/tests/test_metafor_alignment.py
index c8040ed..760f386 100644
--- a/pymare/tests/test_metafor_alignment.py
+++ b/pymare/tests/test_metafor_alignment.py
@@ -101,13 +101,12 @@ def assert_matches(results, case):
 
     Notes
     -----
-    The p-value is compared only under the two corrected options. Under
-    ``"wald"`` PyMARE computes it as ``1 - |0.5 - Phi(z)| * 2``, which cancels
-    catastrophically in the far tail and disagrees with metafor's
-    ``2 * pnorm(-|z|)`` by an arbitrarily large *relative* amount once the p-value
-    drops below about 1e-15. That predates this branch and is not what this module
-    measures; the adjustment's own path goes through ``2 * t.sf(...)``, which does
-    not cancel and agrees to 3e-15.
+    The p-value is compared under all three options. It used to be checked only
+    under the two corrected ones, because ``"wald"`` computed it as
+    ``1 - |0.5 - Phi(z)| * 2``, which cancels catastrophically in the far tail
+    and disagreed with metafor's ``2 * pnorm(-|z|)`` by an arbitrarily large
+    *relative* amount below about 1e-15. Both paths now go through a log
+    survival function, so both agree with metafor at the tolerance above.
     """
     stats = results.get_fe_stats()
     for key, expected in (
@@ -117,8 +116,7 @@ def assert_matches(results, case):
         ("ci_u", "ci_ub"),
     ):
         assert np.allclose(np.ravel(stats[key]), case[expected], rtol=RTOL, atol=ATOL), key
-    if case["test"] != "z":
-        assert np.allclose(np.ravel(stats["p"]), case["pval"], rtol=RTOL, atol=ATOL)
+    assert np.allclose(np.ravel(stats["p"]), case["pval"], rtol=RTOL, atol=ATOL)
 
     if case["dof"] is None:
         assert results.fe_dof is None
diff --git a/pymare/tests/test_results.py b/pymare/tests/test_results.py
index 035ff02..1435a68 100644
--- a/pymare/tests/test_results.py
+++ b/pymare/tests/test_results.py
@@ -4,6 +4,7 @@
 
 import numpy as np
 import pytest
+import scipy.stats as ss
 
 from pymare import Dataset
 from pymare.estimators import (
@@ -133,9 +134,10 @@ def test_mrr_get_fe_stats(results):
     """Test MetaRegressionResults.get_fe_stats."""
     stats = results.get_fe_stats()
     assert isinstance(stats, dict)
-    assert set(stats.keys()) == {"est", "se", "ci_l", "ci_u", "z", "p"}
+    assert set(stats.keys()) == {"est", "se", "ci_l", "ci_u", "z", "p", "logp"}
     assert np.allclose(stats["ci_l"].T, [-7.4651, -1.9693], atol=1e-4)
     assert np.allclose(stats["p"].T, [0.9728, 0.5186], atol=1e-4)
+    assert np.allclose(stats["logp"], np.log(stats["p"]))
     # A t reference with K - P = 6 degrees of freedom, not a normal one.
     assert np.all(results.fe_dof == 6)
 
@@ -160,15 +162,28 @@ def test_mrr_get_heterogeneity_stats(results_2d):
     assert round(stats["I^2"][0], 4) == 88.8487
     assert round(stats["H"][0], 4) == 2.9946
     assert stats["p(Q)"][0] < 1e-5
+    assert np.allclose(stats["logp(Q)"], np.log(stats["p(Q)"]))
 
 
 def test_mrr_to_df(results):
     """Test conversion of MetaRegressionResults to DataFrame."""
     df = results.to_df()
-    assert df.shape == (2, 7)
-    col_names = {"estimate", "p-value", "z-score", "ci_0.025", "ci_0.975", "se", "name"}
+    assert df.shape == (2, 8)
+    col_names = {
+        "estimate",
+        "p-value",
+        "-log10(p)",
+        "z-score",
+        "ci_0.025",
+        "ci_0.975",
+        "se",
+        "name",
+    }
     assert set(df.columns) == col_names
     assert np.allclose(df["p-value"].values, [0.9728, 0.5186], atol=1e-4)
+    # The table reports base 10; get_fe_stats reports natural logs.
+    assert np.allclose(df["-log10(p)"].values, -np.log10(df["p-value"].values))
+    assert np.allclose(results.get_fe_stats()["logp"].ravel(), np.log(df["p-value"].values))
 
 
 def test_small_variance_mrr_to_df(small_variance_results, small_variance_dataset):
@@ -186,10 +201,20 @@ def test_small_variance_mrr_to_df(small_variance_results, small_variance_dataset
     release before the adjustment reported.
     """
     df = small_variance_results.to_df()
-    assert df.shape == (2, 7)
-    col_names = {"estimate", "p-value", "z-score", "ci_0.025", "ci_0.975", "se", "name"}
+    assert df.shape == (2, 8)
+    col_names = {
+        "estimate",
+        "p-value",
+        "-log10(p)",
+        "z-score",
+        "ci_0.025",
+        "ci_0.975",
+        "se",
+        "name",
+    }
     assert set(df.columns) == col_names
     assert np.all(np.isnan(df["p-value"].values))
+    assert np.all(np.isnan(df["-log10(p)"].values))
     assert np.all(small_variance_results.fe_se == 0.0)
 
     unadjusted = (
@@ -197,9 +222,15 @@ def test_small_variance_mrr_to_df(small_variance_results, small_variance_dataset
         .fit_dataset(small_variance_dataset)
         .summary()
     )
-    assert np.allclose(
-        unadjusted.to_df()["p-value"].values, [1, np.finfo(np.float64).eps], atol=1e-4
-    )
+    # 8.2e-23, not the machine epsilon this used to report. Two separate losses
+    # produced that number: ``1 - |0.5 - Phi(z)| * 2`` cancels to exactly 0 by
+    # z = 9.8, and the epsilon floor then presented that 0 as 2.2e-16.
+    unadjusted_df = unadjusted.to_df()
+    assert np.allclose(unadjusted_df["p-value"].values, [1.0, 8.2097e-23], rtol=1e-4)
+    # 22.09 reads as "1e-22" at a glance, which is what the column is for. A
+    # p-value of exactly 1 has to come out at positive zero, not -0.0.
+    assert np.allclose(unadjusted_df["-log10(p)"].values, [0.0, 22.0857], atol=1e-4)
+    assert not np.signbit(unadjusted_df["-log10(p)"].values[0])
 
 
 def test_estimator_summary(dataset):
@@ -352,7 +383,7 @@ def test_heterogeneity_is_undefined_when_the_design_exhausts_the_df():
     results = WeightedLeastSquares().fit_dataset(dataset).summary()
 
     stats = results.get_heterogeneity_stats()
-    assert all(np.all(np.isnan(stats[key])) for key in ("Q", "p(Q)", "I^2", "H"))
+    assert all(np.all(np.isnan(stats[key])) for key in ("Q", "p(Q)", "logp(Q)", "I^2", "H"))
 
 
 def test_undefined_standard_errors_do_not_read_as_significant():
@@ -418,7 +449,7 @@ def test_group_weighted_heterogeneity_matches_collapsed_reference():
         .get_heterogeneity_stats()
     )
 
-    for key in ("Q", "p(Q)", "I^2", "H"):
+    for key in ("Q", "p(Q)", "logp(Q)", "I^2", "H"):
         assert np.allclose(observed[key], expected[key])
 
 
@@ -561,3 +592,117 @@ def test_combination_permutation_freezes_the_correlation_the_estimator_used():
         estimated.permutation_test(n_perm=4).perm_p["fe_p"],
         supplied.permutation_test(n_perm=4).perm_p["fe_p"],
     )
+
+
+# -----------------------------------------------------------------------------
+# The tail, in logs
+# -----------------------------------------------------------------------------
+
+#: The z-score a p-value floored at ``numpy.finfo(float).eps`` maps back to.
+#: Every estimate more significant than 2.2e-16 used to be reported at exactly
+#: this value, whatever its actual evidence.
+EPSILON_CEILING = ss.norm.isf(np.finfo(np.float64).eps / 2)
+
+
+#: Effect sizes spanning the old ceiling: the first two land above p = 2.2e-16,
+#: the last three below it and so used to be reported identically.
+EFFECT_SWEEP = (0.6, 3.0, 100.0, 1e4, 1e40)
+
+
+def test_knapp_hartung_z_climbs_past_the_epsilon_ceiling(extreme_effect_results):
+    """The reported z has to keep separating estimates after p passes 2.2e-16.
+
+    The transform from the t tail to a normal deviate went through a p-value
+    floored at machine epsilon, so it saturated at ``norm.isf(eps / 2) = 8.21``.
+    Estimates whose evidence differed by forty orders of magnitude all came back
+    as 8.21, and a z-map thresholded above that selected nothing.
+    """
+    z = np.array(
+        [np.ravel(extreme_effect_results(e).get_fe_stats()["z"])[0] for e in EFFECT_SWEEP]
+    )
+
+    assert np.all(np.diff(z) > 0)
+    assert (z > EPSILON_CEILING).sum() == 3
+    assert z[-1] > 2 * EPSILON_CEILING
+
+
+def test_z_and_p_describe_the_same_tail(extreme_effect_results):
+    """Consistency between the two is the whole reason z is transformed at all.
+
+    A threshold on z and the corresponding threshold on p must select the same
+    estimates, which makes this a round trip through ``ndtri_exp`` and back. It
+    is checked on the log scale because these p-values are far too small for a
+    relative comparison to survive on the linear one.
+    """
+    for effect in EFFECT_SWEEP:
+        stats = extreme_effect_results(effect).get_fe_stats()
+        assert np.allclose(
+            ss.norm.logsf(np.abs(stats["z"])) + np.log(2), stats["logp"], rtol=1e-12
+        )
+        assert np.allclose(stats["p"], np.exp(stats["logp"]))
+
+
+def test_logp_outlives_p_under_a_normal_reference(extreme_effect_results):
+    """A two-tailed normal p is exactly zero from |z| = 38.5 onwards.
+
+    That is a limit of the representation, not of the evidence, and it is
+    reached on data no more extreme than a well-powered fixed-effects fit.
+    """
+    stats = extreme_effect_results(10.0, correction="wald").get_fe_stats()
+
+    assert np.ravel(stats["p"])[0] == 0.0
+    assert np.allclose(np.ravel(stats["logp"])[0], -9926.1741, rtol=1e-8)
+    # Under a normal reference the statistic needs no transform, so it is
+    # reported as it stands and agrees with metafor's zval exactly.
+    assert np.allclose(stats["z"], stats["est"] / stats["se"])
+
+
+def test_to_df_stays_readable_where_the_p_value_column_is_zero(extreme_effect_results):
+    """The summary table is the only view some callers use, so it needs the log.
+
+    A ``p-value`` of 0.0 beside an ``estimate`` of 10 is the reported bug in its
+    other form: the number is a limit of the column, not of the evidence.
+    """
+    df = extreme_effect_results(10.0, correction="wald").to_df()
+
+    assert df["p-value"].values[0] == 0.0
+    assert np.allclose(df["-log10(p)"].values[0], 4310.8826, rtol=1e-8)
+
+
+def test_heterogeneity_logp_outlives_p_of_q(dataset):
+    """Q grows with the number of observations, so its tail underflows too."""
+    y, v = dataset.y, dataset.y * 0 + 1.0
+    wide = Dataset(y=np.tile(y, (60, 1)), v=np.tile(v, (60, 1)))
+    stats = DerSimonianLaird().fit_dataset(wide).summary().get_heterogeneity_stats()
+
+    assert np.ravel(stats["p(Q)"])[0] == 0.0
+    assert np.allclose(np.ravel(stats["logp(Q)"])[0], -1654.11549, rtol=1e-8)
+
+
+def test_combination_results_rebuild_from_logp_but_not_from_p():
+    """The log is the primitive of the three, and the container has to prefer it.
+
+    A combined z of 60 is unremarkable for a few hundred inputs, and its
+    one-tailed p-value is about 1e-785 -- zero, in a double. Reconstructing from
+    that zero gives an infinite z; reconstructing from the log does not.
+    """
+    z = np.array([[60.0]])
+    logp = ss.norm.logsf(z)
+
+    from_logp = CombinationTestResults(None, None, logp=logp)
+    assert np.allclose(from_logp.z, z)
+    assert np.ravel(from_logp.p)[0] == 0.0
+
+    from_z = CombinationTestResults(None, None, z=z)
+    assert np.allclose(from_z.logp, logp)
+
+    # The lossy direction, asserted so that it stays a documented limit rather
+    # than a surprise: p alone cannot recover what it has already lost.
+    from_p = CombinationTestResults(None, None, p=np.exp(logp))
+    assert not np.isfinite(from_p.z).all()
+
+
+def test_combination_results_still_require_one_array():
+    """Adding a third way to construct one must not make all three optional."""
+    with pytest.raises(ValueError, match="One of 'z', 'p' or 'logp'"):
+        CombinationTestResults(None, None)
diff --git a/pymare/tests/test_stats.py b/pymare/tests/test_stats.py
index a7fd1bd..139f4dc 100644
--- a/pymare/tests/test_stats.py
+++ b/pymare/tests/test_stats.py
@@ -1,9 +1,11 @@
 """Tests for pymare.stats."""
 
 import warnings
+from unittest import mock
 
 import numpy as np
 import pytest
+import scipy.stats as ss
 
 from pymare import stats
 from pymare.estimators import DerSimonianLaird
@@ -24,6 +26,7 @@
     estimate_null_correlation,
     group_mean,
     knapp_hartung_cov_and_dof,
+    log_chi2_sf,
     normalize_group_weights,
     one_sample_t_from_sufficient_statistics,
     satterthwaite_dof,
@@ -1073,3 +1076,77 @@ def test_undo_centering_shrinkage_handles_several_blocks(block_correlation, cent
         block = recovered[start : start + size, start : start + size]
         assert np.allclose(block[~np.eye(size, dtype=bool)], rho, atol=1e-6)
         start += size
+
+
+#: Reference values for :func:`~pymare.stats.log_chi2_sf`, computed with mpmath
+#: at 60 decimal digits as ``log(gammainc(df/2, q/2, inf, regularized=True))``.
+#: Pinned rather than recomputed because mpmath is not a dependency, and taken
+#: from mpmath rather than from SciPy because the point of the function is the
+#: range where SciPy returns ``-inf``. The first three are inside that range and
+#: check the function against a working reference; the rest are outside it.
+LOG_CHI2_SF_REFERENCE = [
+    (3.84, 1.0, -2.994862227180027),
+    (10.0, 4.0, -3.208240530771945),
+    (538.0522161241299, 79.0, -158.26903055960412),
+    (3228.313296744782, 478.0, -924.0499398956518),
+    (5380.522161241301, 798.0, -1535.4587062578848),
+    (1e5, 478.0, -48492.943841385284),
+    (1e6, 2.0, -500000.0),
+]
+
+
+@pytest.mark.parametrize("q,df,expected", LOG_CHI2_SF_REFERENCE)
+def test_log_chi2_sf_matches_arbitrary_precision(q, df, expected):
+    """The tail has to stay accurate long after a double-precision p is zero."""
+    assert np.allclose(log_chi2_sf(q, df), expected, rtol=1e-13)
+
+
+def test_log_chi2_sf_defers_to_scipy_wherever_scipy_can_answer():
+    """Nothing is gained by replacing an answer SciPy already gets right."""
+    df = 7.0
+    q = np.linspace(0.1, 700.0, 2000)
+    assert np.all(np.isfinite(ss.chi2.logsf(q, df)))  # the premise of the sweep
+    assert np.array_equal(log_chi2_sf(q, df), ss.chi2.logsf(q, df))
+
+
+def test_log_chi2_sf_converges_far_inside_its_limit():
+    """The iteration cap has to be generous where the fraction is actually used.
+
+    Convergence slows as ``q`` approaches ``df`` from above -- 722 iterations at
+    ``df = 1e6``, ``q = df + 2`` -- which is why the fraction is confined to the
+    region SciPy cannot reach. This pins that the confinement works: across six
+    orders of magnitude of ``df``, the first ``q`` whose tail underflows already
+    needs single-digit iterations.
+    """
+    for df in (2.0, 478.0, 1e4, 1e6):
+        # The smallest q whose double-precision tail is zero, to a few digits.
+        q = np.array([float(x) for x in np.geomspace(df + 2.0, 1e4 * (df + 2.0), 4000)])
+        underflowed = q[ss.chi2.sf(q, df) == 0.0]
+        assert underflowed.size, df
+
+        edge = underflowed[0]
+        # Converged means moving the cap cannot move the answer.
+        with mock.patch.object(stats, "_GAMMA_CF_MAX_ITER", 8):
+            capped = log_chi2_sf(edge, df)
+        assert np.isfinite(capped)
+        assert np.allclose(capped, log_chi2_sf(edge, df), rtol=1e-14)
+
+
+def test_log_chi2_sf_handles_degenerate_inputs():
+    """An infinite statistic has an exactly zero tail; NaN has no tail at all."""
+    result = log_chi2_sf([np.inf, np.nan, 0.0], 4.0)
+
+    assert result[0] == -np.inf
+    assert np.isnan(result[1])
+    assert result[2] == 0.0
+
+
+def test_log_chi2_sf_broadcasts_and_stays_monotone():
+    """One df per column is how get_heterogeneity_stats calls it."""
+    q = np.array([[10.0, 500.0], [3000.0, 1e5]])
+    df = np.array([4.0, 478.0])
+    logp = log_chi2_sf(q, df)
+
+    assert logp.shape == q.shape
+    assert np.all(np.diff(logp, axis=0) < 0)
+    assert np.allclose(logp[0, 0], log_chi2_sf(10.0, 4.0))