diff --git a/src/decijax/__init__.py b/src/decijax/__init__.py
index b0393d6..872fcac 100644
--- a/src/decijax/__init__.py
+++ b/src/decijax/__init__.py
@@ -4,6 +4,7 @@
     acquisition_functions,
     acquisition_maximizer,
     decision_maker,
+    maths,
     models,
     search_space,
     test_functions,
@@ -15,6 +16,7 @@
     "acquisition_functions",
     "acquisition_maximizer",
     "decision_maker",
+    "maths",
     "models",
     "search_space",
     "test_functions",
diff --git a/src/decijax/acquisition_functions/__init__.py b/src/decijax/acquisition_functions/__init__.py
index 1886a70..e7e8358 100644
--- a/src/decijax/acquisition_functions/__init__.py
+++ b/src/decijax/acquisition_functions/__init__.py
@@ -8,6 +8,7 @@
 )
 from decijax.acquisition_functions.expected_improvement import (
     ExpectedImprovement,
+    LogExpectedImprovement,
 )
 from decijax.acquisition_functions.probability_of_improvement import (
     LogProbabilityOfImprovement,
@@ -20,6 +21,7 @@
     "AbstractAcquisitionFunctionBuilder",
     "AbstractSinglePointAcquisitionFunctionBuilder",
     "ExpectedImprovement",
+    "LogExpectedImprovement",
     "LogProbabilityOfImprovement",
     "SinglePointAcquisitionFunction",
     "ThompsonSampling",
diff --git a/src/decijax/acquisition_functions/expected_improvement.py b/src/decijax/acquisition_functions/expected_improvement.py
index 1c4dd1f..9e55eef 100644
--- a/src/decijax/acquisition_functions/expected_improvement.py
+++ b/src/decijax/acquisition_functions/expected_improvement.py
@@ -3,6 +3,7 @@
 from collections.abc import Mapping
 
 import jax.numpy as jnp
+from jax.nn import logmeanexp
 from jax.scipy.stats import norm
 from jaxtyping import (
     Array,
@@ -13,6 +14,10 @@
     AbstractSinglePointAcquisitionFunctionBuilder,
     SinglePointAcquisitionFunction,
 )
+from decijax.maths import (
+    _log1mexp,
+    _log_abs_z_cdf_div_pdf,
+)
 from decijax.models import (
     ProbabilisticModel,
     SupportsGaussianPrediction,
@@ -24,12 +29,73 @@
 )
 
 
+def _log_ei_helper(z: Float[Array, "..."]) -> Float[Array, "..."]:
+    r"""Compute $`\log(\phi(z) + z \Phi(z))`$ in a numerically stable manner.
+
+    Evaluated directly this underflows to zero, gradient included, once
+    $`z \lesssim -6`$ in double precision. Implements the piecewise formulation of
+    [Ament et al., 2023](https://arxiv.org/abs/2310.20708), Eq. 9:
+
+    ```math
+    \begin{cases}
+    \log\left(\phi(z) + z \Phi(z)\right) & z > -1 \\
+    -\frac{z^2}{2} - c_1 + \text{log1mexp}\left(
+        \log\left(\text{erfcx}\left(\frac{-z}{\sqrt{2}}\right)|z|\right) + c_2
+    \right) & \tau < z \leq -1 \\
+    -\frac{z^2}{2} - c_1 - 2\log|z| & z \leq \tau
+    \end{cases}
+    ```
+
+    where $`c_1 = \frac{\log(2\pi)}{2}`$ and $`c_2 = \frac{\log(\pi / 2)}{2}`$. The
+    threshold $`\tau`$ between the last two branches is $`-10^6`$ in double
+    precision and $`-10^3`$ in single, following BoTorch. Each branch clamps on
+    *both* sides of its boundaries to keep gradients clean.
+
+    Args:
+        z: Scaled improvement.
+
+    Returns:
+        $`\log(\phi(z) + z \Phi(z))`$, elementwise.
+    """
+    if z.dtype == jnp.float64:
+        eps_bound = -1e6
+    elif z.dtype == jnp.float32:
+        eps_bound = -1e3
+    else:
+        raise NotImplementedError(
+            f"LogExpectedImprovement does not support dtype {z.dtype}."
+        )
+
+    log_sqrt_2pi = 0.5 * jnp.log(2.0 * jnp.pi)
+
+    # Branch 1 (z > -1): direct computation.
+    z_upper = jnp.maximum(z, -1.0)
+    direct = jnp.log(norm.pdf(z_upper) + z_upper * norm.cdf(z_upper))
+
+    # Branch 2 (eps_bound < z <= -1): stable computation using erfcx.
+    z_mid = jnp.clip(z, eps_bound, -1.0)
+    w = _log_abs_z_cdf_div_pdf(z_mid)
+    stable = -0.5 * z_mid**2 - log_sqrt_2pi + _log1mexp(w)
+
+    # Branch 3 (z <= eps_bound): asymptotic, where `_log1mexp` can no longer resolve w.
+    z_lower = jnp.minimum(z, eps_bound)
+    asymptotic = -0.5 * z_lower**2 - log_sqrt_2pi - 2.0 * jnp.log(-z_lower)
+
+    return jnp.where(
+        z > -1.0,
+        direct,
+        jnp.where(z > eps_bound, stable, asymptotic),
+    )
+
+
 class ExpectedImprovement(AbstractSinglePointAcquisitionFunctionBuilder):
     """Standard Expected Improvement acquisition function.
 
-    As introduced by [Močkus, 1974](https://link.springer.com/chapter/10.1007/3-540-07165-2_55). The "best" incumbent value is defined as the highest posterior
-    mean value evaluated at the previously observed points. This enables the
-    acquisition function to be utilised with noisy observations.
+    As introduced by
+    [Močkus, 1974](https://link.springer.com/chapter/10.1007/3-540-07165-2_55). The
+    "best" incumbent value is defined as the highest posterior mean value evaluated
+    at the previously observed points. This enables the acquisition function to be
+    utilised with noisy observations.
     """
 
     def build_acquisition_function(
@@ -41,9 +107,12 @@ def build_acquisition_function(
 
         This computes the expected improvement over the "best" of the previously
         observed points, utilising the posterior distribution of the surrogate model.
-        For posterior distribution $`f(\cdot)`$, and best incumbent value $`\eta`$, this is defined as:
+        For posterior distribution $`f(\cdot)`$, and best incumbent value $`\eta`$,
+        this is defined as:
+
         ```math
-        \alpha_{\text{EI}}(\mathbf{x}) = \mathbb{E}\left[\max(0, f(\mathbf{x}) - \eta)\right]
+        \alpha_{\text{EI}}(\mathbf{x})
+        = \mathbb{E}\left[\max(0, f(\mathbf{x}) - \eta)\right]
         ```
 
         For models carrying a leading sample axis (e.g. fully Bayesian GPs), the
@@ -85,3 +154,82 @@ def _expected_improvement(x: Float[Array, "N D"]) -> Float[Array, "N 1"]:
             return jnp.mean(ei, axis=0)[:, None]  # marginalise over S -> [N, 1]
 
         return _expected_improvement
+
+
+class LogExpectedImprovement(AbstractSinglePointAcquisitionFunctionBuilder):
+    r"""Numerically stable Log Expected Improvement acquisition function [1].
+
+    Given a predictive posterior distribution of the objective function $f$, the log
+    expected improvement at a test point $x$ is defined as:
+
+    $$\text{LogEI}(x) = \log \mathbb{E}\left[\max(0, f(x) - f(x^*))\right]$$
+
+    where $x^*$ is the maximiser of the posterior mean at previously observed values
+    (to handle noisy observations).
+
+    Being a strictly increasing transform of the expected improvement, this shares
+    its maximiser exactly, but is far better behaved as an optimisation target:
+    expected improvement vanishes to *exactly* zero, gradient included, once the
+    scaled improvement falls below roughly $-40$ in double precision, and those flat
+    regions come to dominate the search space in higher dimensions.
+
+    References:
+    ----------
+    [1] Ament, S., Daulton, S., Eriksson, D., Balandat, M., & Bakshy, E. (2023).
+    Unexpected improvement to expected improvement for Bayesian optimization.
+    Advances in Neural Information Processing Systems, 36.
+    """
+
+    def build_acquisition_function(
+        self,
+        models: Mapping[str, ProbabilisticModel],
+        key: KeyArray,
+    ) -> SinglePointAcquisitionFunction:
+        r"""Build the Log Expected Improvement acquisition function.
+
+        The expected improvement factorises as $`\sigma \cdot h(z)`$, for scaled
+        improvement $`z = \frac{\mu - \eta}{\sigma}`$ and
+        $`h(z) = \phi(z) + z\Phi(z)`$, so that:
+
+        ```math
+        \alpha_{\text{LogEI}}(\mathbf{x}) = \log \sigma(\mathbf{x}) + \log h(z)
+        ```
+
+        with the second term computed by `_log_ei_helper`. For models carrying a
+        leading sample axis (e.g. fully Bayesian GPs), it is computed per sample and
+        reduced with a log-mean-exp, the correct marginalisation
+        $`\log \mathbb{E}_\theta[\alpha_{\text{EI},\theta}(\mathbf{x})]`$.
+
+        Args:
+            models: Dictionary of models used to form the acquisition function. One
+                model must correspond to the `OBJECTIVE` key and support Gaussian
+                prediction, as we use the objective posterior to form the acquisition
+                function.
+            key: JAX PRNG key used for random number generation. Since the log
+                expected improvement is computed deterministically, the key is not
+                used.
+
+        Returns:
+            The Log Expected Improvement acquisition function to be *maximised* in
+            order to decide which point to query next.
+        """
+        self.check_objective_present(models)
+        objective_model = models[OBJECTIVE]
+
+        if not isinstance(objective_model, SupportsGaussianPrediction):
+            raise ValueError(
+                "Objective model must support Gaussian prediction to compute the "
+                "Log Expected Improvement."
+            )
+
+        eta = get_best_latent_observation_val(objective_model)  # [S, 1]
+
+        def _log_expected_improvement(x: Float[Array, "N D"]) -> Float[Array, "N 1"]:
+            latent_dist = objective_model.predict(x)
+            mean = latent_dist.mean  # [S, N]
+            std = latent_dist.stddev  # [S, N]
+            z = (mean - eta) / std
+            log_ei = _log_ei_helper(z) + jnp.log(std)  # [S, N]
+            return logmeanexp(log_ei, axis=0)[:, None]  # marginalise over S -> [N, 1]
+
+        return _log_expected_improvement
diff --git a/src/decijax/maths.py b/src/decijax/maths.py
new file mode 100644
index 0000000..d72e019
--- /dev/null
+++ b/src/decijax/maths.py
@@ -0,0 +1,100 @@
+"""A collection of numerically stable mathematical functions."""
+
+import jax.numpy as jnp
+from jax.scipy.special import erfc
+from jaxtyping import (
+    Array,
+    Float,
+)
+
+
+def _log1mexp(x: Float[Array, "..."]) -> Float[Array, "..."]:
+    r"""Compute $`\log(1 - e^x)`$ for $`x < 0`$ without catastrophic cancellation.
+
+    Which form is accurate depends on how close $`e^x`$ is to one, hence the split
+    at $`-\log 2`$. See [Mächler,
+    2012](https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf).
+
+    Args:
+        x: Strictly negative array.
+
+    Returns:
+        $`\log(1 - e^x)`$, elementwise.
+    """
+    log_2 = jnp.log(2.0)
+    # Each branch is clamped to its own side of the split, so the discarded one
+    # cannot poison the gradient of the surviving one.
+    return jnp.where(
+        x > -log_2,
+        jnp.log(-jnp.expm1(jnp.maximum(x, -log_2))),
+        jnp.log1p(-jnp.exp(jnp.minimum(x, -log_2))),
+    )
+
+
+def _erfcx(x: Float[Array, "..."]) -> Float[Array, "..."]:
+    r"""Compute the scaled complementary error function $`e^{x^2}\text{erfc}(x)`$.
+
+    We need to define our own erfcx function, because `jax.scipy.special.erfcx` has
+    a bug where it currently returns zero on $`[26.54, 26.64]`$ (double) and
+    $`[9.19, 9.42]`$ (single). Once
+    [jax-ml/jax#38607](https://github.com/jax-ml/jax/issues/38607) is fixed, replace
+    this function with the native JAX implementation.
+
+    Args:
+        x: Strictly positive array.
+
+    Returns:
+        $`e^{x^2}\text{erfc}(x)`$, elementwise.
+    """
+    # Beyond this point the asymptotic expansion below has reached machine precision.
+    # It must stay *below* the point at which erfc underflows (26.54 in double, 9.19
+    # in single), which is exactly the bug JAX's own implementation has.
+    if x.dtype == jnp.float64:
+        asymptotic_bound = 13.7
+    elif x.dtype == jnp.float32:
+        asymptotic_bound = 6.9
+    else:
+        raise NotImplementedError(f"_erfcx does not support dtype {x.dtype}.")
+
+    # erfcx ~ (sqrt(pi) x)^-1 sum_k c_k x^-2k, with c_k = (-1)^k (2k-1)!! / 2^k, in
+    # descending degree for `jnp.polyval`.
+    coeffs = jnp.asarray(
+        [
+            7918.06640625,
+            -1055.7421875,
+            162.421875,
+            -29.53125,
+            6.5625,
+            -1.875,
+            0.75,
+            -0.5,
+            1.0,
+        ],
+        dtype=x.dtype,
+    )
+
+    large = x > asymptotic_bound
+    # Clamped so the discarded branch can neither overflow nor poison the gradient.
+    safe_x = jnp.where(large, jnp.ones_like(x), x)
+    direct = jnp.exp(jnp.square(safe_x)) * erfc(safe_x)
+    asymptotic = jnp.polyval(coeffs, 1.0 / jnp.square(x)) / (x * jnp.sqrt(jnp.pi))
+    return jnp.where(large, asymptotic, direct)
+
+
+def _log_abs_z_cdf_div_pdf(z: Float[Array, "..."]) -> Float[Array, "..."]:
+    r"""Compute $`\log(|z| \Phi(z) / \phi(z))`$ for $`z < 0`$.
+
+    Deep in the tail $`\Phi(z)`$ and $`\phi(z)`$ both underflow to zero even though
+    their ratio is perfectly ordinary: at $`z = -40`$ both are around $`10^{-349}`$,
+    while $`|z| \Phi(z) / \phi(z)`$ is $`0.999`$. `_erfcx` cancels the shared
+    $`e^{-z^2/2}`$ factor analytically, so neither tiny number is ever formed.
+
+    Args:
+        z: Strictly negative array.
+
+    Returns:
+        Elementwise; strictly negative, approaching zero from below.
+    """
+    neg_inv_sqrt_2 = -(2.0**-0.5)
+    log_sqrt_pi_div_2 = 0.5 * jnp.log(jnp.pi / 2.0)
+    return jnp.log(_erfcx(neg_inv_sqrt_2 * z) * jnp.abs(z)) + log_sqrt_pi_div_2
diff --git a/tests/test_acquisition_functions/test_acquisition_functions.py b/tests/test_acquisition_functions/test_acquisition_functions.py
index 9019bab..2c751f8 100644
--- a/tests/test_acquisition_functions/test_acquisition_functions.py
+++ b/tests/test_acquisition_functions/test_acquisition_functions.py
@@ -5,6 +5,7 @@
 )
 from decijax.acquisition_functions.expected_improvement import (
     ExpectedImprovement,
+    LogExpectedImprovement,
 )
 from decijax.acquisition_functions.probability_of_improvement import (
     ProbabilityOfImprovement,
@@ -29,7 +30,12 @@
 )  # Sampling with tfp causes JAX to raise a UserWarning due to some internal logic around jnp.argsort
 @pytest.mark.parametrize(
     "acquisition_function_builder",
-    [ExpectedImprovement, ProbabilityOfImprovement, ThompsonSampling],
+    [
+        ExpectedImprovement,
+        LogExpectedImprovement,
+        ProbabilityOfImprovement,
+        ThompsonSampling,
+    ],
 )
 def test_acquisition_function_no_objective_model_raises_error(
     acquisition_function_builder: type[AbstractSinglePointAcquisitionFunctionBuilder],
@@ -49,7 +55,12 @@ def test_acquisition_function_no_objective_model_raises_error(
 )  # Sampling with tfp causes JAX to raise a UserWarning due to some internal logic around jnp.argsort
 @pytest.mark.parametrize(
     "acquisition_function_builder",
-    [ExpectedImprovement, ProbabilityOfImprovement, ThompsonSampling],
+    [
+        ExpectedImprovement,
+        LogExpectedImprovement,
+        ProbabilityOfImprovement,
+        ThompsonSampling,
+    ],
 )
 def test_model_without_required_capability_raises_error(
     acquisition_function_builder: type[AbstractSinglePointAcquisitionFunctionBuilder],
@@ -66,7 +77,12 @@ def test_model_without_required_capability_raises_error(
 
 @pytest.mark.parametrize(
     "acquisition_function_builder",
-    [ExpectedImprovement, ProbabilityOfImprovement, ThompsonSampling],
+    [
+        ExpectedImprovement,
+        LogExpectedImprovement,
+        ProbabilityOfImprovement,
+        ThompsonSampling,
+    ],
 )
 @pytest.mark.parametrize(
     "test_target_function",
diff --git a/tests/test_acquisition_functions/test_expected_improvement.py b/tests/test_acquisition_functions/test_expected_improvement.py
index 2760413..480ce03 100644
--- a/tests/test_acquisition_functions/test_expected_improvement.py
+++ b/tests/test_acquisition_functions/test_expected_improvement.py
@@ -1,9 +1,12 @@
+import jax
 import jax.numpy as jnp
 import jax.random as jr
 import numpyro.distributions as dist
 import pytest
 from decijax.acquisition_functions.expected_improvement import (
     ExpectedImprovement,
+    LogExpectedImprovement,
+    _log_ei_helper,
 )
 from decijax.models import GPJaxConjugateGP
 from decijax.test_functions.continuous_functions import (
@@ -48,3 +51,145 @@ def test_expected_improvement_acquisition_function_correct_values(
     mc_ei = jnp.expand_dims(jnp.mean(jnp.maximum(samples - eta, 0), 0), -1)
     assert jnp.all(ei >= 0)
     assert jnp.allclose(ei, mc_ei, rtol=0.03, atol=1e-6)
+
+
+@pytest.mark.parametrize(
+    "test_target_function",
+    [NegativeForrester(), NegativeLogarithmicGoldsteinPrice()],
+)
+@pytest.mark.parametrize("key", [jr.key(42), jr.key(10)])
+def test_log_expected_improvement_acquisition_function_correct_values(
+    test_target_function: AbstractContinuousTestFunction,
+    key: KeyArray,
+):
+    # LogEI must be the log of the (marginalised) EI: exp(LogEI) should recover both
+    # the analytic EI and its Monte-Carlo estimate.
+    data_key, ei_acq_key, log_ei_acq_key, test_key, mc_key = jr.split(key, 5)
+    dataset = test_target_function.generate_dataset(num_points=10, key=data_key)
+    posterior = generate_dummy_conjugate_posterior(dataset, test_target_function)
+    model = GPJaxConjugateGP(posterior=posterior, dataset=dataset)
+    models = {OBJECTIVE: model}
+    ei_fn = ExpectedImprovement().build_acquisition_function(models, ei_acq_key)
+    log_ei_fn = LogExpectedImprovement().build_acquisition_function(
+        models, log_ei_acq_key
+    )
+    test_x = test_target_function.generate_test_points(100, test_key)
+    ei = ei_fn(test_x)
+    log_ei = log_ei_fn(test_x)
+    latent_dist = posterior.predict(test_x, dataset)
+    latent_mean = latent_dist.mean
+    latent_var = latent_dist.variance
+    samples = dist.Normal(loc=latent_mean, scale=jnp.sqrt(latent_var)).sample(
+        mc_key, sample_shape=(10000,)
+    )
+    eta = get_best_latent_observation_val(model)
+    mc_ei = jnp.expand_dims(jnp.mean(jnp.maximum(samples - eta, 0), 0), -1)
+    assert log_ei.shape == (100, 1)
+    assert jnp.all(jnp.isfinite(log_ei))
+    assert jnp.allclose(jnp.exp(log_ei), ei, rtol=1e-6, atol=1e-6)
+    assert jnp.allclose(jnp.exp(log_ei), mc_ei, rtol=0.03, atol=1e-6)
+
+
+def test_log_ei_helper_is_stable_in_the_tails():
+    # The naive log(phi(z) + z * Phi(z)) underflows to -inf, gradient included, well
+    # within the range a maximiser explores. Straddles both branch points, -1 and -1e6.
+    z = jnp.array(
+        [
+            -1e100,
+            -1e10,
+            -1e6 - 1.0,
+            -1e6 + 1.0,
+            -1e3,
+            -300.0,
+            -100.0,
+            -37.6,  # see test_log_ei_helper_is_accurate_where_erfcx_is_broken
+            -30.0,
+            -10.0,
+            -3.0,
+            -1.0 - 1e-6,
+            -1.0 + 1e-6,
+            0.0,
+            1.0,
+            10.0,
+        ]
+    )
+    log_ei = _log_ei_helper(z)
+    grad = jax.grad(lambda z_: jnp.sum(_log_ei_helper(z_)))(z)
+
+    assert jnp.all(jnp.isfinite(log_ei))
+    assert jnp.all(jnp.isfinite(grad))
+    # LogEI is strictly increasing in the scaled improvement.
+    assert jnp.all(grad > 0)
+    assert jnp.all(jnp.diff(log_ei) > 0)
+
+    # Where the naive form still has headroom in float64, the two must agree.
+    naive_z = z[z >= -30.0]
+    naive = jnp.log(
+        jax.scipy.stats.norm.pdf(naive_z) + naive_z * jax.scipy.stats.norm.cdf(naive_z)
+    )
+    assert jnp.allclose(_log_ei_helper(naive_z), naive, rtol=1e-6, atol=1e-6)
+
+    # And the naive form is genuinely unusable where the helper is not.
+    deep_tail = jnp.array([-100.0, -1e3])
+    assert jnp.all(
+        jnp.log(
+            jax.scipy.stats.norm.pdf(deep_tail)
+            + deep_tail * jax.scipy.stats.norm.cdf(deep_tail)
+        )
+        == -jnp.inf
+    )
+
+
+def test_log_ei_helper_matches_asymptotic_expansion_in_the_tails():
+    # Finite and monotone is not enough; the deep tail must be *correct*, and there
+    # the naive form has underflowed and is no use as a reference. For large |z|,
+    # h(z) = phi(z) * (1/z^2 - 3/z^4 + 15/z^6 - ...).
+    z = -jnp.concatenate(
+        [
+            jnp.linspace(10.0, 1000.0, 400),
+            jnp.array([1e4, 1e6, 1e10, 1e50]),  # straddles the -1e6 branch point
+        ]
+    )
+    y = 1.0 / z**2
+    series = 1.0 - 3.0 * y + 15.0 * y**2 - 105.0 * y**3 + 945.0 * y**4
+    expected = -0.5 * z**2 - 0.5 * jnp.log(2.0 * jnp.pi) + jnp.log(y * series)
+
+    # Compared in absolute terms, against values ranging down to -1e100.
+    assert jnp.allclose(_log_ei_helper(z), expected, rtol=0.0, atol=1e-4)
+
+
+def test_log_ei_helper_is_accurate_where_erfcx_is_broken():
+    # The band where jax.scipy.special.erfcx returns 0.0 maps to z in
+    # [-37.68, -37.54]. Going through it naively gives a silently wrong (finite,
+    # monotone) value, so pin it down against a reworking that reintroduces the bug.
+    z = jnp.linspace(-37.7, -37.5, 501)
+    y = 1.0 / z**2
+    series = 1.0 - 3.0 * y + 15.0 * y**2 - 105.0 * y**3 + 945.0 * y**4
+    expected = -0.5 * z**2 - 0.5 * jnp.log(2.0 * jnp.pi) + jnp.log(y * series)
+
+    assert jnp.allclose(_log_ei_helper(z), expected, rtol=0.0, atol=1e-4)
+
+
+@pytest.mark.parametrize(
+    "test_target_function",
+    [NegativeForrester(), NegativeLogarithmicGoldsteinPrice()],
+)
+def test_log_expected_improvement_is_jit_and_grad_safe(
+    test_target_function: AbstractContinuousTestFunction,
+):
+    # The maximiser jits and differentiates the closure, so it must be pure.
+    key = jr.key(42)
+    data_key, acq_key, test_key = jr.split(key, 3)
+    dataset = test_target_function.generate_dataset(num_points=10, key=data_key)
+    posterior = generate_dummy_conjugate_posterior(dataset, test_target_function)
+    model = GPJaxConjugateGP(posterior=posterior, dataset=dataset)
+    log_ei_fn = LogExpectedImprovement().build_acquisition_function(
+        {OBJECTIVE: model}, acq_key
+    )
+    test_x = test_target_function.generate_test_points(10, test_key)
+
+    jitted = jax.jit(log_ei_fn)(test_x)
+    grad = jax.grad(lambda x: jnp.sum(log_ei_fn(x)))(test_x)
+
+    assert jnp.allclose(jitted, log_ei_fn(test_x))
+    assert jnp.all(jnp.isfinite(grad))
diff --git a/tests/test_maths.py b/tests/test_maths.py
new file mode 100644
index 0000000..53eccb3
--- /dev/null
+++ b/tests/test_maths.py
@@ -0,0 +1,42 @@
+import jax
+import jax.numpy as jnp
+import numpy as np
+import scipy.special
+from decijax.maths import _erfcx
+
+
+def test_erfcx_matches_scipy_across_its_whole_range():
+    # We own `_erfcx`, so it needs its own coverage. These straddle the asymptotic
+    # branch point (13.7 in double), the band JAX gets wrong, and the extremes where
+    # the direct form would overflow or underflow.
+    x = jnp.array(
+        [
+            1e-3,
+            1.0,
+            6.9,
+            13.0,
+            13.7 - 1e-6,
+            13.7 + 1e-6,
+            20.0,
+            26.6,
+            30.0,
+            100.0,
+            1e4,
+            1e7,
+        ]
+    )
+    assert jnp.allclose(_erfcx(x), scipy.special.erfcx(np.asarray(x)), rtol=1e-13)
+    assert jnp.all(_erfcx(x) > 0.0)
+    assert jnp.all(jnp.isfinite(jax.grad(lambda x_: jnp.sum(_erfcx(x_)))(x)))
+
+    # Dense, since sparse points step straight over a narrow band of wrong values.
+    x_dense = jnp.linspace(1e-3, 40.0, 20001)
+    assert jnp.allclose(
+        _erfcx(x_dense), scipy.special.erfcx(np.asarray(x_dense)), rtol=1e-13
+    )
+
+
+def test_jax_erfcx_is_still_broken_so_the_vendored_one_is_still_needed():
+    # If this ever fails, jax-ml/jax#38607 is fixed and `_erfcx` can go.
+    x = jnp.linspace(26.55, 26.64, 101)
+    assert jnp.all(jax.scipy.special.erfcx(x) == 0.0)