diff --git a/docs/api-reference/index.md b/docs/api-reference/index.md
index 196442178..deb3d33f9 100644
--- a/docs/api-reference/index.md
+++ b/docs/api-reference/index.md
@@ -65,5 +65,6 @@ and possible confusion of `theta` (from Bragg’s law) with `theta` in spherical
    logging
    metadata
    peaks
+   smoothing
    tof
 ```
diff --git a/src/scippneutron/__init__.py b/src/scippneutron/__init__.py
index 4f883682c..717f64b6d 100644
--- a/src/scippneutron/__init__.py
+++ b/src/scippneutron/__init__.py
@@ -32,6 +32,7 @@
     'data',
     'metadata',
     'peaks',
+    'smoothing',
     'tof',
 ]
 
diff --git a/src/scippneutron/smoothing.py b/src/scippneutron/smoothing.py
new file mode 100644
index 000000000..ad1440ed6
--- /dev/null
+++ b/src/scippneutron/smoothing.py
@@ -0,0 +1,621 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2026 Scipp contributors (https://github.com/scipp)
+"""Kernel smoothing of one-dimensional sampled signals.
+
+``smooth`` applies a translation-invariant displacement distribution, for a kernel
+whose width is constant in coordinate units. ``smooth_relative`` instead applies a
+distribution of fractional displacements, so its width scales with the coordinate.
+
+Both procedures accept irregularly spaced samples. They interpolate onto a uniform
+or geometric working grid that is no coarser than the tightest input spacing,
+integrate the kernel probability over grid cells, renormalize the convolution at
+finite boundaries, and interpolate the result back to the original coordinates.
+"""
+
+from __future__ import annotations
+
+from numbers import Real
+from typing import Any, Protocol, TypeAlias, cast
+
+import numpy as np
+import scipp as sc
+from numpy.typing import ArrayLike, NDArray
+from scipy.signal import convolve
+from scipy.stats import norm, triang, uniform
+
+__all__ = [
+    "smooth",
+    "smooth_relative",
+]
+
+
+class _Distribution(Protocol):
+    def cdf(self, x: ArrayLike) -> Any: ...
+
+    def ppf(self, probability: ArrayLike) -> Any: ...
+
+    def support(self) -> tuple[Any, Any]: ...
+
+
+_Kernel: TypeAlias = str | _Distribution
+_FloatArray: TypeAlias = NDArray[np.float64]
+_IntArray: TypeAlias = NDArray[np.int64]
+_ScippArray: TypeAlias = sc.DataArray | sc.Variable
+
+_BUILTIN_KERNELS: dict[str, _Distribution] = {
+    # Standard Gaussian.
+    "gaussian": norm(),
+    "normal": norm(),
+    # Centered rectangle on [-1, 1].
+    "rectangle": uniform(loc=-1.0, scale=2.0),
+    "rect": uniform(loc=-1.0, scale=2.0),
+    "box": uniform(loc=-1.0, scale=2.0),
+    "boxcar": uniform(loc=-1.0, scale=2.0),
+    "uniform": uniform(loc=-1.0, scale=2.0),
+    # Centered triangle on [-1, 1], peak at 0.
+    # Users can pass triang(c=...) themselves for asymmetric triangles.
+    "triangle": triang(c=0.5, loc=-1.0, scale=2.0),
+    "triangular": triang(c=0.5, loc=-1.0, scale=2.0),
+}
+
+
+def _as_kernel_distribution(kernel: _Kernel) -> _Distribution:
+    if isinstance(kernel, str):
+        try:
+            return _BUILTIN_KERNELS[kernel.lower()]
+        except KeyError:
+            valid = ", ".join(sorted(_BUILTIN_KERNELS))
+            raise ValueError(
+                f"unknown kernel {kernel!r}; expected one of: {valid}"
+            ) from None
+
+    required = ("cdf", "ppf", "support")
+    missing = [name for name in required if not hasattr(kernel, name)]
+    if missing:
+        raise TypeError(
+            "kernel must be a scipy.stats distribution-like object with methods "
+            f"{', '.join(required)}; missing {', '.join(missing)}"
+        )
+
+    return kernel
+
+
+def _kernel_support(dist: _Distribution) -> tuple[float, float]:
+    try:
+        support_min, support_max = dist.support()
+    except TypeError as e:
+        raise TypeError(
+            "kernel must be a fully specified distribution. For distributions "
+            "with shape parameters, pass a frozen distribution such as "
+            "triang(c=0.5), not triang."
+        ) from e
+    return float(support_min), float(support_max)
+
+
+def _trim_kernel_weights(
+    offsets: _IntArray, weights: _FloatArray
+) -> tuple[_IntArray, _FloatArray]:
+    nonzero = np.flatnonzero(weights > 0.0)
+    if nonzero.size == 0:
+        return offsets, weights
+
+    # Trim zero-only ends, but preserve offset zero for convolution alignment.
+    zero = -offsets[0]
+    first = min(nonzero[0], zero)
+    last = max(nonzero[-1], zero) + 1
+    weights = weights[first:last]
+    return offsets[first:last], weights / weights.sum()
+
+
+class _Geometry(Protocol):
+    """
+    A grid on which the kernel is translation invariant.
+
+    Both smoothing procedures displace a coordinate x by a distribution Z. The
+    displacement law differs, but in both cases there is a working coordinate
+    ``u`` in which the displacement is an additive offset independent of x, so
+    that a single set of weights applies at every grid point. ``offset`` and
+    ``displacement`` convert between a displacement z and its offset in u; they
+    are inverses of each other.
+    """
+
+    #: Names the grid in user-facing messages.
+    name: str
+
+    def check(self, x: _FloatArray) -> None:
+        """Reject coordinates outside the domain of the displacement law."""
+
+    def coordinate(self, x: _FloatArray) -> _FloatArray:
+        """Working coordinate u(x)."""
+
+    def points(self, start: float, stop: float, count: int) -> _FloatArray:
+        """``count`` samples from ``start`` to ``stop``, uniform in u."""
+
+    def displacement_min(self, scale: float) -> float:
+        """Smallest displacement z that keeps x within the valid domain."""
+
+    def offset(self, scale: float, z: ArrayLike) -> Any:
+        """Offset in u produced by displacement z."""
+
+    def displacement(self, scale: float, u: ArrayLike) -> Any:
+        """Displacement z producing an offset u."""
+
+
+class _Uniform:
+    """Constant kernel width: ``x' = x + scale * Z``, additive in x itself."""
+
+    name = "uniform"
+
+    def check(self, x: _FloatArray) -> None:
+        pass
+
+    def coordinate(self, x: _FloatArray) -> _FloatArray:
+        return x
+
+    def points(self, start: float, stop: float, count: int) -> _FloatArray:
+        return np.linspace(start, stop, count)
+
+    def displacement_min(self, scale: float) -> float:
+        return -np.inf
+
+    def offset(self, scale: float, z: ArrayLike) -> Any:
+        return scale * np.asarray(z, dtype=float)
+
+    def displacement(self, scale: float, u: ArrayLike) -> Any:
+        return np.asarray(u, dtype=float) / scale
+
+
+class _Geometric:
+    """
+    Relative kernel width: ``x' = x * (1 + scale * Z)``, additive in ``log(x)``.
+
+    Equivalently, the kernel is
+
+        K(x, x') = 1 / (scale * x) * f((x' - x) / (scale * x))
+
+    for a distribution with PDF f. Displacements are restricted to ``x' > 0``,
+    that is ``1 + scale * z > 0``.
+    """
+
+    name = "geometric"
+
+    def check(self, x: _FloatArray) -> None:
+        if np.any(x <= 0):
+            raise ValueError("x must be positive")
+
+    def coordinate(self, x: _FloatArray) -> _FloatArray:
+        return cast(_FloatArray, np.log(x))
+
+    def points(self, start: float, stop: float, count: int) -> _FloatArray:
+        return np.geomspace(start, stop, count)
+
+    def displacement_min(self, scale: float) -> float:
+        return -1.0 / scale
+
+    def offset(self, scale: float, z: ArrayLike) -> Any:
+        # log1p(-1) is -inf, which _kernel_weights treats as unbounded support.
+        with np.errstate(divide="ignore", invalid="ignore"):
+            return np.log1p(scale * np.asarray(z, dtype=float))
+
+    def displacement(self, scale: float, u: ArrayLike) -> Any:
+        return np.expm1(np.asarray(u, dtype=float)) / scale
+
+
+def _kernel_weights(
+    geometry: _Geometry,
+    spacing: float,
+    scale: float,
+    kernel: _Kernel,
+    tail: float,
+    max_offset: int,
+) -> tuple[_IntArray, _FloatArray]:
+    """
+    Cell-integrated kernel weights on a grid of the given spacing in u.
+
+    Weights are exact integrals of the kernel probability over grid cells, so
+    the result is a proper quadrature of the smoothing integral rather than a
+    point sampling of the kernel.
+    """
+    dist = _as_kernel_distribution(kernel)
+
+    z_min = geometry.displacement_min(scale)
+    p_min = float(dist.cdf(z_min))
+    reachable_mass = 1.0 - p_min
+    if not np.isfinite(reachable_mass) or reachable_mass <= 0.0:
+        raise ValueError("kernel has no mass in the valid domain for this scale")
+
+    support_min, support_max = _kernel_support(dist)
+    u_left = float(geometry.offset(scale, max(z_min, support_min)))
+    u_right = float(geometry.offset(scale, support_max))
+
+    # Unbounded support, or support reaching the edge of the valid domain,
+    # gives an infinite bound in u. Truncate at the requested tail instead.
+    if not (np.isfinite(u_left) and np.isfinite(u_right)):
+        probabilities = (
+            p_min + np.array([0.5 * tail, 1.0 - 0.5 * tail]) * reachable_mass
+        )
+        u_left, u_right = (
+            float(value) for value in geometry.offset(scale, dist.ppf(probabilities))
+        )
+        if not np.isfinite(u_right):
+            raise ValueError("right kernel bound is not finite; increase tail")
+
+    # Cells are centered at m*h and span [(m-1/2)h, (m+1/2)h].
+    # Include offset zero even for one-sided kernels and clamp the stencil to
+    # offsets that can contribute to the finite input.
+    m_min = int(np.clip(np.floor(u_left / spacing + 0.5), -max_offset, 0))
+    m_max = int(np.clip(np.ceil(u_right / spacing - 0.5), 0, max_offset))
+    m = np.arange(m_min, m_max + 1, dtype=np.int64)
+
+    lower = geometry.displacement(scale, (m - 0.5) * spacing)
+    upper = geometry.displacement(scale, (m + 0.5) * spacing)
+    weights = np.maximum(dist.cdf(upper) - dist.cdf(lower), 0.0)
+    return _trim_kernel_weights(m, weights)
+
+
+def _valid_weight_sums(n: int, m: _IntArray, w: _FloatArray) -> _FloatArray:
+    """
+    Boundary normalization.
+
+    Equivalent to
+
+        np.convolve(np.ones(n), w[::-1], mode="full")[start:start+n]
+
+    but O(n), not another full convolution.
+    """
+    i = np.arange(n, dtype=np.int64)
+
+    m_min = int(m[0])
+    m_max = int(m[-1])
+
+    lower = np.maximum(m_min, -i)
+    upper = np.minimum(m_max, n - 1 - i)
+
+    cumsum = np.empty(w.size + 1, dtype=float)
+    cumsum[0] = 0.0
+    np.cumsum(w, out=cumsum[1:])
+
+    return cast(
+        _FloatArray,
+        cumsum[upper - m_min + 1] - cumsum[lower - m_min],
+    )
+
+
+def _smooth_with_weights(
+    y: _FloatArray, offsets: _IntArray, weights: _FloatArray
+) -> _FloatArray:
+    # Desired operation:
+    #
+    #   out[i] = sum_m weights[m] * y[i + m]
+    #
+    # scipy convolution reverses the second argument, hence weights[::-1].
+    # FFT convolution would spread a single NaN or infinity over the entire
+    # output. SciPy recommends the direct method for non-finite inputs.
+    method = "auto" if np.all(np.isfinite(y)) else "direct"
+    full = cast(_FloatArray, convolve(y, weights[::-1], mode="full", method=method))
+    start = int(offsets[-1])
+    numerator = full[start : start + y.size]
+    denominator = _valid_weight_sums(y.size, offsets, weights)
+
+    out = np.full_like(numerator, np.nan)
+    np.divide(numerator, denominator, out=out, where=denominator > 0.0)
+    return out
+
+
+def _validate_max_grid_points(max_grid_points: int) -> None:
+    if isinstance(max_grid_points, bool | np.bool_) or not isinstance(
+        max_grid_points, int | np.integer
+    ):
+        raise TypeError("max_grid_points must be an integer")
+    if max_grid_points < 2:
+        raise ValueError("max_grid_points must be at least 2")
+
+
+def _smooth_values(
+    geometry: _Geometry,
+    x: ArrayLike,
+    y: ArrayLike,
+    scale: float,
+    kernel: _Kernel,
+    tail: float,
+    max_grid_points: int,
+) -> _FloatArray:
+    x = np.asarray(x, dtype=float)
+    y = np.asarray(y, dtype=float)
+
+    if x.ndim != 1 or y.ndim != 1:
+        raise ValueError("x and y must be one-dimensional")
+    if x.size != y.size:
+        raise ValueError("x and y must have the same length")
+    if not np.isfinite(scale) or scale < 0:
+        raise ValueError("scale must be non-negative")
+    if not np.isfinite(tail) or not (0.0 < tail < 1.0):
+        raise ValueError("tail must be between 0 and 1")
+    _validate_max_grid_points(max_grid_points)
+    if np.any(~np.isfinite(x)):
+        raise ValueError("x must contain only finite values")
+    geometry.check(x)
+    if np.any(np.diff(x) <= 0):
+        raise ValueError("x must be strictly increasing")
+    if x.size < 2 or scale == 0:
+        return y.copy()
+
+    u = geometry.coordinate(x)
+    du = np.diff(u)
+    u_range = float(u[-1] - u[0])
+
+    # Preserve an existing regular grid. Otherwise choose a grid at least as
+    # dense as the smallest input spacing, measured in the working coordinate.
+    if np.allclose(du, du[0], rtol=1e-7, atol=0.0):
+        k = float(x.size)
+    else:
+        # Coordinates that are distinct but collide, or nearly so, in the
+        # working coordinate give an infinite point count. Let it propagate to
+        # the guard below rather than raising here.
+        with np.errstate(divide="ignore", over="ignore"):
+            k = np.ceil(u_range / np.min(du)) + 1.0
+
+    if not (k <= max_grid_points) and k > 2 * x.size:
+        raise ValueError(
+            f"{geometry.name} resampling would require too many points, exceeding "
+            f"max_grid_points={max_grid_points:,}. Increase max_grid_points to "
+            "allow a larger grid."
+        )
+    k = int(k)
+
+    spacing = u_range / (k - 1)
+    if not np.isfinite(spacing) or spacing <= 0:
+        raise ValueError("grid spacing must be positive")
+
+    offsets, weights = _kernel_weights(
+        geometry,
+        spacing=spacing,
+        scale=scale,
+        kernel=kernel,
+        tail=tail,
+        max_offset=k - 1,
+    )
+    xp = geometry.points(float(x[0]), float(x[-1]), k)
+    zg = _smooth_with_weights(np.interp(xp, x, y), offsets, weights)
+    return cast(_FloatArray, np.interp(x, xp, zg))
+
+
+def _scipp_input(
+    x: object, y: object
+) -> tuple[sc.Variable, sc.Variable, sc.DataArray | None]:
+    """Validate untyped user input and reduce it to a coordinate/data pair."""
+    template: sc.DataArray | None = None
+    if isinstance(x, sc.DataArray):
+        if y is not None:
+            raise TypeError("y must be omitted when x is a DataArray")
+        if x.ndim != 1:
+            raise sc.DimensionError("data must be one-dimensional")
+        if x.dim not in x.coords:
+            raise sc.CoordError("data must have a dimension coordinate")
+        if x.coords.is_edges(x.dim):
+            raise sc.CoordError("the dimension coordinate must not contain bin edges")
+        if x.masks:
+            raise ValueError("smoothing data with masks is not supported")
+        template = x
+        y = x.data
+        x = x.coords[x.dim]
+    elif not isinstance(x, sc.Variable) or not isinstance(y, sc.Variable):
+        raise TypeError("expected a DataArray or a pair of Variables")
+
+    if x.ndim != 1 or y.ndim != 1:
+        raise sc.DimensionError("x and y must be one-dimensional")
+    if x.dims != y.dims:
+        raise sc.DimensionError("x and y must have the same dimension")
+    if x.is_binned or y.is_binned:
+        raise sc.DTypeError("x and y must not be binned")
+    if y.variances is not None:
+        raise sc.VariancesError(
+            "Smoothing signals with variances is not supported because it would "
+            "introduce correlations between data points."
+        )
+    return x, y, template
+
+
+def _scipp_output(
+    template: sc.DataArray | None, y: sc.Variable, values: _FloatArray
+) -> _ScippArray:
+    data = sc.array(dims=y.dims, values=values, unit=y.unit)
+    if template is None:
+        return data
+
+    # The new container shares unchanged coordinates, but its data is independent.
+    out = template.copy(deep=False)
+    out.data = data
+    return out
+
+
+def _scale_in_coordinate_unit(scale: sc.Variable, x: sc.Variable) -> float:
+    if not isinstance(scale, sc.Variable):
+        raise TypeError("scale must be a scipp.Variable")
+    if scale.ndim != 0:
+        raise sc.DimensionError("scale must be a scalar")
+    if scale.variances is not None:
+        raise sc.VariancesError("kernel scales with variances are not supported")
+    return float(scale.to(unit=x.unit).value)
+
+
+def _dimensionless_scale(scale: object) -> float:
+    if not isinstance(scale, sc.Variable):
+        if not isinstance(scale, Real):
+            raise TypeError("scale must be a real number or a scipp.Variable")
+        return float(scale)
+    if scale.ndim != 0:
+        raise sc.DimensionError("scale must be a scalar")
+    if scale.variances is not None:
+        raise sc.VariancesError("kernel scales with variances are not supported")
+    return float(scale.to(unit=sc.units.dimensionless).value)
+
+
+def smooth(
+    x: _ScippArray,
+    y: sc.Variable | None = None,
+    *,
+    scale: sc.Variable,
+    kernel: _Kernel = "gaussian",
+    tail: float = 1e-12,
+    max_grid_points: int = 1_000_000,
+) -> _ScippArray:
+    """Smooth sampled data with a translation-invariant kernel.
+
+    The kernel describes a distribution of displacements ``Z``, with displaced
+    coordinates given by ``x' = x + scale * Z``. At the boundaries, the kernel
+    is renormalized over the available finite input domain.
+
+    Input that is not uniformly spaced is interpolated to a uniform grid,
+    smoothed, and interpolated back to the original coordinates.
+
+    Parameters
+    ----------
+    x:
+        One-dimensional data to smooth, or strictly increasing sample
+        coordinates. A data array must have a dimension coordinate.
+    y:
+        Values to smooth when ``x`` contains the sample coordinates. Must be a
+        one-dimensional variable with the same dimension as ``x``. Must be
+        omitted when ``x`` is a data array.
+    scale:
+        Scale factor for the displacement distribution. Must be a scalar with a
+        unit compatible with the coordinate. Set to zero to return a copy of the
+        input without smoothing.
+    kernel:
+        Kernel distribution. The canonical names are ``'gaussian'``,
+        ``'boxcar'``, and ``'triangular'``. They represent a standard normal
+        distribution, a uniform distribution on [-1, 1], and a symmetric
+        triangular distribution on [-1, 1], respectively. Other aliases are
+        accepted. Alternatively, provide a fully specified distribution with
+        ``cdf``, ``ppf``, and ``support`` methods.
+    tail:
+        Total probability omitted when truncating a kernel with unbounded
+        support.
+    max_grid_points:
+        Intermediate uniform grids no larger than this are always allowed. Larger
+        grids may be rejected to guard against excessive resampling.
+
+    Returns
+    -------
+    :
+        Smoothed data of the same type as the input. Coordinates and units are
+        preserved.
+
+    Raises
+    ------
+    ValueError
+        If the inputs have invalid values, if a data array has masks, if a string
+        does not identify a supported kernel, or if the required intermediate
+        grid exceeds ``max_grid_points``.
+    scipp.DimensionError
+        If the inputs are not one-dimensional, a pair of variables does not
+        have matching dimensions, or ``scale`` is not scalar.
+    scipp.CoordError
+        If a data array has no dimension coordinate or has a bin-edge
+        coordinate.
+    scipp.UnitError
+        If the unit of ``scale`` is incompatible with the coordinate unit.
+    scipp.VariancesError
+        If the signal or ``scale`` has variances.
+    TypeError
+        If ``scale`` is not a variable, ``kernel`` is not a distribution-like
+        object, or ``max_grid_points`` is not an integer.
+    """
+    x, y, template = _scipp_input(x, y)
+    values = _smooth_values(
+        _Uniform(),
+        x.values,
+        y.values,
+        scale=_scale_in_coordinate_unit(scale, x),
+        kernel=kernel,
+        tail=tail,
+        max_grid_points=max_grid_points,
+    )
+    return _scipp_output(template, y, values)
+
+
+def smooth_relative(
+    x: _ScippArray,
+    y: sc.Variable | None = None,
+    *,
+    scale: float | sc.Variable,
+    kernel: _Kernel = "gaussian",
+    tail: float = 1e-12,
+    max_grid_points: int = 1_000_000,
+) -> _ScippArray:
+    """Smooth sampled data with a kernel of relative width.
+
+    The kernel describes a distribution of relative displacements ``Z``, with
+    displaced coordinates given by ``x' = x * (1 + scale * Z)``. At the
+    boundaries, the kernel is renormalized over the available finite input
+    domain.
+
+    Input that is not geometrically spaced is interpolated to a geometric grid,
+    smoothed, and interpolated back to the original coordinates.
+
+    Parameters
+    ----------
+    x:
+        One-dimensional data to smooth, or positive, strictly increasing
+        one-dimensional sample coordinates. A data array must have a
+        dimension coordinate.
+    y:
+        Values to smooth when ``x`` contains the sample coordinates. Must be a
+        one-dimensional variable with the same dimension as ``x``. Must be
+        omitted when ``x`` is a data array.
+    scale:
+        Dimensionless scale factor for the relative-displacement distribution.
+        May be a real number or a scalar, dimensionless variable. Set to zero to
+        return a copy of the input without smoothing.
+    kernel:
+        Kernel distribution. The canonical names are ``'gaussian'``,
+        ``'boxcar'``, and ``'triangular'``. They represent a standard normal
+        distribution, a uniform distribution on [-1, 1], and a symmetric
+        triangular distribution on [-1, 1], respectively. Other aliases are
+        accepted. Alternatively, provide a fully specified distribution with
+        ``cdf``, ``ppf``, and ``support`` methods.
+    tail:
+        Total probability omitted when truncating a kernel with unbounded support
+        or support reaching the nonpositive coordinate domain.
+    max_grid_points:
+        Intermediate geometric grids no larger than this are always allowed. Larger
+        grids may be rejected to guard against excessive resampling.
+
+    Returns
+    -------
+    :
+        Smoothed data of the same type as the input. Coordinates and units are
+        preserved.
+
+    Raises
+    ------
+    ValueError
+        If the inputs have invalid values, if a data array has masks, if a
+        string does not identify a supported kernel, or if the required
+        intermediate grid exceeds ``max_grid_points``.
+    scipp.DimensionError
+        If the inputs are not one-dimensional or a pair of variables does not
+        have matching dimensions, or ``scale`` is not scalar.
+    scipp.CoordError
+        If a data array has no dimension coordinate or has a bin-edge
+        coordinate.
+    scipp.UnitError
+        If ``scale`` is a variable with a non-dimensionless unit.
+    scipp.VariancesError
+        If the signal or ``scale`` has variances.
+    TypeError
+        If ``scale`` is neither a real number nor a variable, ``kernel`` is not
+        a distribution-like object, or ``max_grid_points`` is not an integer.
+    """
+    x, y, template = _scipp_input(x, y)
+    values = _smooth_values(
+        _Geometric(),
+        x.values,
+        y.values,
+        scale=_dimensionless_scale(scale),
+        kernel=kernel,
+        tail=tail,
+        max_grid_points=max_grid_points,
+    )
+    return _scipp_output(template, y, values)
diff --git a/tests/smoothing_test.py b/tests/smoothing_test.py
new file mode 100644
index 000000000..0c5f8f5c0
--- /dev/null
+++ b/tests/smoothing_test.py
@@ -0,0 +1,836 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2026 Scipp contributors (https://github.com/scipp)
+
+import numpy as np
+import pytest
+import scipp as sc
+from scipy.special import ndtr
+from scipy.stats import uniform
+
+import scippneutron as scn
+from scippneutron.smoothing import smooth, smooth_relative
+
+
+def _normal_pdf(z):
+    z = np.asarray(z, dtype=float)
+    return np.exp(-0.5 * z**2) / np.sqrt(2.0 * np.pi)
+
+
+def _quadratic(x):
+    return 1.0 + 0.3 * x + 0.7 * x**2
+
+
+def _exact_smoothed_quadratic_with_sigma(x, sigma, lower, upper):
+    x = np.asarray(x, dtype=float)
+    sigma = np.asarray(sigma, dtype=float)
+
+    z_lower = (lower - x) / sigma
+    z_upper = (upper - x) / sigma
+    normalization = ndtr(z_upper) - ndtr(z_lower)
+
+    mean_z = (_normal_pdf(z_lower) - _normal_pdf(z_upper)) / normalization
+    mean_z_squared = (
+        1.0
+        + (z_lower * _normal_pdf(z_lower) - z_upper * _normal_pdf(z_upper))
+        / normalization
+    )
+
+    mean_x = x + sigma * mean_z
+    mean_x_squared = x**2 + 2.0 * x * sigma * mean_z + sigma**2 * mean_z_squared
+    return 1.0 + 0.3 * mean_x + 0.7 * mean_x_squared
+
+
+def _exact_smoothed_quadratic(x, alpha, lower, upper):
+    """Exactly smooth ``_quadratic`` with a relative finite-domain Gaussian."""
+    return _exact_smoothed_quadratic_with_sigma(x, alpha * np.asarray(x), lower, upper)
+
+
+def _geometric_cell_centers(lower, upper, size):
+    log_step = np.log(upper / lower) / size
+    return lower * np.exp((np.arange(size) + 0.5) * log_step)
+
+
+def _linear_cell_centers(lower, upper, size):
+    step = (upper - lower) / size
+    return lower + (np.arange(size) + 0.5) * step
+
+
+def _variables(x, y):
+    return (
+        sc.array(dims=['x'], values=np.asarray(x), unit='m'),
+        sc.array(
+            dims=['x'],
+            values=np.asarray(y),
+            unit='counts',
+        ),
+    )
+
+
+def _smooth_values(smooth, x, y, **kwargs):
+    x, y = _variables(x, y)
+    return smooth(x, y, **kwargs).values
+
+
+def _quadratic_smoothing_error(*, size, alpha, tail=1e-9):
+    lower = 0.1
+    upper = 0.9
+    x = _geometric_cell_centers(lower, upper, size)
+    actual = _smooth_values(
+        smooth_relative,
+        x,
+        _quadratic(x),
+        scale=alpha,
+        tail=tail,
+    )
+    expected = _exact_smoothed_quadratic(x, alpha, lower, upper)
+    return actual - expected
+
+
+def _fixed_width_quadratic_smoothing_error(*, size, width, tail=1e-9):
+    lower = -0.4
+    upper = 0.9
+    x = _linear_cell_centers(lower, upper, size)
+    x_var, y_var = _variables(x, _quadratic(x))
+    actual = smooth(
+        x_var,
+        y_var,
+        scale=sc.scalar(width, unit='m'),
+        tail=tail,
+    ).values
+    expected = _exact_smoothed_quadratic_with_sigma(x, width, lower, upper)
+    return actual - expected
+
+
+def test_fixed_width_gaussian_matches_exact_quadratic_on_regular_grid():
+    error = _fixed_width_quadratic_smoothing_error(size=4000, width=0.1)
+
+    assert np.max(np.abs(error)) < 5e-7
+    assert np.sqrt(np.mean(error**2)) < 1e-7
+
+
+def test_fixed_width_gaussian_error_is_second_order_in_grid_spacing():
+    sizes = np.array([500, 1000, 2000, 4000])
+    errors = np.array(
+        [
+            np.max(np.abs(_fixed_width_quadratic_smoothing_error(size=size, width=0.1)))
+            for size in sizes
+        ]
+    )
+
+    observed_orders = np.log2(errors[:-1] / errors[1:])
+    np.testing.assert_allclose(observed_orders, 2.0, atol=0.06)
+
+
+def test_fixed_width_gaussian_error_scales_with_tail_until_grid_error_dominates():
+    tails = np.array([1e-2, 1e-3, 1e-4, 1e-5])
+    errors = np.array(
+        [
+            np.max(
+                np.abs(
+                    _fixed_width_quadratic_smoothing_error(
+                        size=4000, width=0.1, tail=tail
+                    )
+                )
+            )
+            for tail in tails
+        ]
+    )
+
+    reduction_per_decade = errors[:-1] / errors[1:]
+    assert np.all((5.0 < reduction_per_decade) & (reduction_per_decade < 15.0))
+
+    grid_limited_errors = np.array(
+        [
+            np.max(
+                np.abs(
+                    _fixed_width_quadratic_smoothing_error(
+                        size=4000, width=0.1, tail=tail
+                    )
+                )
+            )
+            for tail in (1e-9, 1e-12)
+        ]
+    )
+    np.testing.assert_allclose(
+        grid_limited_errors[0], grid_limited_errors[1], rtol=0.05
+    )
+
+
+@pytest.mark.parametrize(
+    ("kernel", "kernel_variance", "max_error"),
+    [
+        ("boxcar", 1.0 / 3.0, 3e-7),
+        ("triangular", 1.0 / 6.0, 5e-8),
+    ],
+)
+def test_fixed_width_compact_kernel_matches_interior_quadratic_moments(
+    kernel, kernel_variance, max_error
+):
+    lower = -0.4
+    upper = 0.9
+    size = 4000
+    width = 0.1
+    x = _linear_cell_centers(lower, upper, size)
+    x_var, y_var = _variables(x, _quadratic(x))
+
+    actual = smooth(
+        x_var,
+        y_var,
+        scale=sc.scalar(width, unit='m'),
+        kernel=kernel,
+    ).values
+    expected = 1.0 + 0.3 * x + 0.7 * (x**2 + width**2 * kernel_variance)
+    interior = (x - width >= lower) & (x + width <= upper)
+
+    assert np.max(np.abs(actual[interior] - expected[interior])) < max_error
+
+
+def test_fixed_width_smoothing_converts_scale_to_coordinate_unit():
+    x = sc.linspace('x', -0.4, 0.9, 100, unit='m')
+    y = sc.array(dims=['x'], values=_quadratic(x.values), unit='counts')
+
+    in_meters = smooth(x, y, scale=sc.scalar(0.1, unit='m'))
+    in_centimeters = smooth(x, y, scale=sc.scalar(10.0, unit='cm'))
+
+    assert sc.identical(in_meters, in_centimeters)
+
+
+def test_fixed_width_smoothing_accepts_distribution():
+    x = sc.linspace('x', -0.4, 0.9, 100, unit='m')
+    y = sc.array(dims=['x'], values=_quadratic(x.values), unit='counts')
+    scale = sc.scalar(0.1, unit='m')
+
+    actual = smooth(
+        x,
+        y,
+        scale=scale,
+        kernel=uniform(loc=-1.0, scale=2.0),
+    )
+    expected = smooth(x, y, scale=scale, kernel="boxcar")
+
+    assert sc.identical(actual, expected)
+
+
+def test_fixed_width_asymmetric_kernel_has_correct_direction():
+    x = sc.arange('x', 0.0, 2.0, 0.1, unit='m')
+    y = sc.array(dims=['x'], values=x.values, unit='counts')
+
+    actual = smooth(
+        x,
+        y,
+        scale=sc.scalar(0.2, unit='m'),
+        kernel=uniform(loc=0.5, scale=1.0),
+    )
+
+    np.testing.assert_allclose(actual.values[3:-3], y.values[3:-3] + 0.2)
+
+
+def test_nonfinite_value_only_affects_overlapping_kernel_windows():
+    size = 10_000
+    x = sc.arange('x', float(size))
+    values = np.ones(size)
+    values[size // 2] = np.nan
+    y = sc.array(dims=['x'], values=values)
+
+    actual = smooth(x, y, scale=sc.scalar(150.0))
+
+    assert np.isfinite(actual.values[0])
+    assert np.isnan(actual.values[size // 2])
+    assert np.isfinite(actual.values[-1])
+
+
+def test_fixed_width_smoothing_accepts_data_array():
+    x = sc.linspace('x', -0.4, 0.9, 100, unit='m')
+    data = sc.DataArray(
+        sc.array(dims=['x'], values=_quadratic(x.values), unit='counts'),
+        coords={'x': x, 'aux': sc.arange('x', 100)},
+    )
+
+    actual = smooth(data, scale=sc.scalar(0.1, unit='m'))
+    expected = smooth(x, data.data, scale=sc.scalar(0.1, unit='m'))
+
+    assert sc.identical(actual.data, expected)
+    assert sc.identical(actual.coords['x'], data.coords['x'])
+    assert sc.identical(actual.coords['aux'], data.coords['aux'])
+
+
+def test_fixed_width_smoothing_resamples_nonuniform_grid():
+    x_values = np.array([0.0, 0.25, 0.75, 1.0])
+    y_values = _quadratic(x_values)
+    scale = sc.scalar(0.2, unit='m')
+    x, y = _variables(x_values, y_values)
+
+    actual = smooth(x, y, scale=scale)
+
+    xp_values = np.linspace(0.0, 1.0, 5)
+    xp, yp = _variables(xp_values, np.interp(xp_values, x_values, y_values))
+    smoothed_yp = smooth(xp, yp, scale=scale)
+    expected = np.interp(x_values, xp_values, smoothed_yp.values)
+    np.testing.assert_allclose(actual.values, expected)
+
+
+def test_fixed_width_gaussian_matches_exact_quadratic_on_nonuniform_grid():
+    lower = -0.4
+    upper = 0.9
+    size = 2000
+    width = 0.1
+    spacing = (upper - lower) / size
+    x = _linear_cell_centers(lower, upper, size)
+    x += 0.2 * spacing * np.sin(np.linspace(0.0, 8.0 * np.pi, size))
+    x_var, y_var = _variables(x, _quadratic(x))
+
+    actual = smooth(
+        x_var,
+        y_var,
+        scale=sc.scalar(width, unit='m'),
+        tail=1e-9,
+    ).values
+
+    intervals = int(np.ceil((x[-1] - x[0]) / np.min(np.diff(x))))
+    resampled_spacing = (x[-1] - x[0]) / intervals
+    expected = _exact_smoothed_quadratic_with_sigma(
+        x,
+        width,
+        x[0] - 0.5 * resampled_spacing,
+        x[-1] + 0.5 * resampled_spacing,
+    )
+
+    assert np.max(np.abs(actual - expected)) < 5e-7
+
+
+def test_fixed_width_smoothing_rejects_resampled_grid_larger_than_limit():
+    x, y = _variables([0.0, 0.125, 1.0], [1.0, 2.0, 3.0])
+
+    with pytest.raises(
+        ValueError,
+        match=r"uniform resampling would require too many points.*max_grid_points=4",
+    ):
+        smooth(x, y, scale=sc.scalar(0.1, unit='m'), max_grid_points=4)
+
+
+def test_fixed_width_smoothing_accepts_resampled_grid_equal_to_limit():
+    x, y = _variables([0.0, 0.25, 1.0], [1.0, 1.0, 1.0])
+
+    actual = smooth(
+        x,
+        y,
+        scale=sc.scalar(0.1, unit='m'),
+        max_grid_points=5,
+    )
+
+    np.testing.assert_allclose(actual.values, y.values)
+
+
+@pytest.mark.parametrize(
+    ("function", "x", "kwargs"),
+    [
+        (
+            smooth,
+            np.linspace(0.0, 1.0, 11),
+            {"scale": sc.scalar(0.1, unit='m')},
+        ),
+        (smooth_relative, np.geomspace(1.0, 2.0, 11), {"scale": 0.1}),
+    ],
+)
+def test_input_grid_larger_than_limit_is_accepted(function, x, kwargs):
+    y = np.ones_like(x)
+
+    actual = _smooth_values(function, x, y, max_grid_points=10, **kwargs)
+
+    np.testing.assert_allclose(actual, y)
+
+
+@pytest.mark.parametrize(
+    ("function", "x", "kwargs"),
+    [
+        (smooth, np.array([0.0, 0.26, 1.0]), {"scale": sc.scalar(0.1, unit='m')}),
+        (smooth_relative, 2.0 ** np.array([0.0, 0.26, 1.0]), {"scale": 0.1}),
+    ],
+)
+def test_modest_resampling_larger_than_limit_is_accepted(function, x, kwargs):
+    y = np.ones_like(x)
+
+    actual = _smooth_values(function, x, y, max_grid_points=4, **kwargs)
+
+    np.testing.assert_allclose(actual, y)
+
+
+@pytest.mark.parametrize(
+    "coordinates",
+    [
+        [1.0, np.nextafter(1.0, 2.0), 2.0],
+        # The spacing ratio overflows to infinity rather than merely being large.
+        [0.0, 5e-324, 1.0],
+        [1.0, np.nextafter(1.0, 2.0), 1e300],
+    ],
+)
+def test_fixed_width_pathologically_close_coordinates_fail_before_allocation(
+    coordinates,
+):
+    x, y = _variables(coordinates, [1.0, 1.0, 1.0])
+
+    with pytest.raises(ValueError, match="exceeding max_grid_points=1,000,000"):
+        smooth(x, y, scale=sc.scalar(0.1, unit='m'))
+
+
+@pytest.mark.parametrize("scale", [np.nan, np.inf, -np.inf, -1.0])
+def test_fixed_width_smoothing_rejects_invalid_scale(scale):
+    x, y = _variables([], [])
+
+    with pytest.raises(ValueError, match="scale must be non-negative"):
+        smooth(x, y, scale=sc.scalar(scale, unit='m'))
+
+
+def test_fixed_width_smoothing_rejects_non_scalar_scale():
+    x, y = _variables([0.0, 1.0], [1.0, 2.0])
+
+    with pytest.raises(sc.DimensionError, match="scale must be a scalar"):
+        smooth(x, y, scale=sc.array(dims=['scale'], values=[0.1], unit='m'))
+
+
+def test_fixed_width_smoothing_rejects_non_variable_scale():
+    x, y = _variables([0.0, 1.0], [1.0, 2.0])
+
+    with pytest.raises(TypeError, match=r"scale must be a scipp\.Variable"):
+        smooth(x, y, scale=0.1)  # type: ignore[arg-type]
+
+
+def test_fixed_width_smoothing_rejects_incompatible_scale_unit():
+    x, y = _variables([0.0, 1.0], [1.0, 2.0])
+
+    with pytest.raises(sc.UnitError):
+        smooth(x, y, scale=sc.scalar(0.1, unit='s'))
+
+
+def test_fixed_width_smoothing_rejects_scale_with_variance():
+    x, y = _variables([0.0, 1.0], [1.0, 2.0])
+
+    with pytest.raises(sc.VariancesError, match="scales with variances"):
+        smooth(
+            x,
+            y,
+            scale=sc.scalar(0.1, variance=0.01, unit='m'),
+        )
+
+
+def test_gaussian_smoothing_matches_exact_quadratic_on_geometric_grid():
+    # The outer cell edges of this grid coincide with the finite integration
+    # domain used by _quadratic_smoothing_error.
+    error = _quadratic_smoothing_error(size=4000, alpha=0.1)
+
+    assert np.max(np.abs(error)) < 5e-7
+    assert np.sqrt(np.mean(error**2)) < 1e-7
+
+
+def test_gaussian_smoothing_error_is_second_order_in_log_grid_spacing():
+    sizes = np.array([500, 1000, 2000, 4000])
+    errors = np.array(
+        [
+            np.max(np.abs(_quadratic_smoothing_error(size=size, alpha=0.1)))
+            for size in sizes
+        ]
+    )
+
+    # Halving the log-grid spacing should reduce midpoint quadrature error by
+    # four, corresponding to second-order convergence.
+    observed_orders = np.log2(errors[:-1] / errors[1:])
+    np.testing.assert_allclose(observed_orders, 2.0, atol=0.06)
+
+
+def test_gaussian_smoothing_error_scales_with_inverse_kernel_width():
+    size = 4000
+    alphas = np.array([0.025, 0.05, 0.1])
+    errors = np.array(
+        [
+            np.max(np.abs(_quadratic_smoothing_error(size=size, alpha=alpha)))
+            for alpha in alphas
+        ]
+    )
+
+    # For a well-resolved narrow Gaussian, the largest error is at a truncated
+    # boundary and scales as h**2 / alpha.
+    log_step = np.log(0.9 / 0.1) / size
+    scaled_errors = errors * alphas / log_step**2
+    assert np.all(np.diff(errors) < 0.0)
+    np.testing.assert_allclose(
+        scaled_errors,
+        np.mean(scaled_errors),
+        rtol=0.08,
+    )
+
+
+def test_gaussian_smoothing_error_scales_with_tail_until_grid_error_dominates():
+    tails = np.array([1e-2, 1e-3, 1e-4, 1e-5])
+    truncation_errors = np.array(
+        [
+            np.max(np.abs(_quadratic_smoothing_error(size=4000, alpha=0.1, tail=tail)))
+            for tail in tails
+        ]
+    )
+
+    # Gaussian moments in the omitted tails add logarithmic factors, so each
+    # decade should improve the result by approximately, but not exactly, ten.
+    reduction_per_decade = truncation_errors[:-1] / truncation_errors[1:]
+    assert np.all((5.0 < reduction_per_decade) & (reduction_per_decade < 15.0))
+
+    # Once tail truncation is negligible, reducing it further cannot improve
+    # the fixed-grid quadrature error.
+    grid_limited_errors = np.array(
+        [
+            np.max(np.abs(_quadratic_smoothing_error(size=4000, alpha=0.1, tail=tail)))
+            for tail in (1e-9, 1e-12)
+        ]
+    )
+    np.testing.assert_allclose(
+        grid_limited_errors[0], grid_limited_errors[1], rtol=0.05
+    )
+
+
+@pytest.mark.parametrize(
+    ("kernel", "relative_variance", "max_error"),
+    [
+        ("boxcar", 1.0 / 3.0, 3e-7),
+        ("triangular", 1.0 / 6.0, 5e-8),
+    ],
+)
+def test_compact_symmetric_kernel_matches_exact_interior_quadratic_moments(
+    kernel, relative_variance, max_error
+):
+    lower = 0.1
+    upper = 0.9
+    size = 4000
+    alpha = 0.1
+    x = _geometric_cell_centers(lower, upper, size)
+
+    actual = _smooth_values(
+        smooth_relative,
+        x,
+        _quadratic(x),
+        scale=alpha,
+        kernel=kernel,
+    )
+    expected = 1.0 + 0.3 * x + 0.7 * x**2 * (1.0 + alpha**2 * relative_variance)
+    interior = (x * (1.0 - alpha) >= lower) & (x * (1.0 + alpha) <= upper)
+
+    assert np.max(np.abs(actual[interior] - expected[interior])) < max_error
+
+
+def test_smoothing_module_is_exposed_by_package():
+    assert scn.smoothing.smooth_relative is smooth_relative
+
+
+def test_accepts_data_array_and_preserves_metadata():
+    x = sc.geomspace('x', 0.1, 0.9, 100, unit='m')
+    data = sc.DataArray(
+        sc.array(dims=['x'], values=_quadratic(x.values), unit='counts'),
+        coords={
+            'x': x,
+            'aux': sc.arange('x', 100, unit='s'),
+            'scalar': sc.scalar(1.2, unit='K'),
+        },
+    )
+    original = data.copy()
+
+    actual = smooth_relative(data, scale=0.1)
+    expected = smooth_relative(x, data.data, scale=0.1)
+
+    assert sc.identical(data, original)
+    assert isinstance(actual, sc.DataArray)
+    assert sc.identical(actual.data, expected)
+    assert sc.identical(actual.coords['x'], data.coords['x'])
+    assert sc.identical(actual.coords['aux'], data.coords['aux'])
+    assert sc.identical(actual.coords['scalar'], data.coords['scalar'])
+
+    actual.values[0] = -1.0
+    assert sc.identical(data, original)
+
+
+def test_relative_smoothing_accepts_dimensionless_variable_scale():
+    x, y = _variables([1.0, 2.0, 3.0], [1.0, 4.0, 9.0])
+
+    actual = smooth_relative(x, y, scale=sc.scalar(0.1))
+    expected = smooth_relative(x, y, scale=0.1)
+
+    assert sc.identical(actual, expected)
+
+
+def test_relative_smoothing_rejects_non_scalar_scale():
+    x, y = _variables([1.0, 2.0], [1.0, 2.0])
+
+    with pytest.raises(sc.DimensionError, match="scale must be a scalar"):
+        smooth_relative(x, y, scale=sc.array(dims=['scale'], values=[0.1]))
+
+
+def test_relative_smoothing_rejects_scale_with_unit():
+    x, y = _variables([1.0, 2.0], [1.0, 2.0])
+
+    with pytest.raises(sc.UnitError):
+        smooth_relative(x, y, scale=sc.scalar(0.1, unit='m'))
+
+
+def test_relative_smoothing_rejects_scale_with_variance():
+    x, y = _variables([1.0, 2.0], [1.0, 2.0])
+
+    with pytest.raises(sc.VariancesError, match="scales with variances"):
+        smooth_relative(x, y, scale=sc.scalar(0.1, variance=0.01))
+
+
+def test_rejects_variable_with_variances():
+    x = sc.geomspace('x', 0.1, 0.9, 100, unit='m')
+    values = _quadratic(x.values)
+    y = sc.array(
+        dims=['x'],
+        values=values,
+        variances=2.0 + x.values,
+        unit='counts',
+    )
+
+    with pytest.raises(sc.VariancesError, match="signals with variances"):
+        smooth_relative(x, y, scale=0.1)
+
+
+def test_rejects_data_array_with_variances():
+    x = sc.geomspace('x', 0.1, 0.9, 100, unit='m')
+    data = sc.DataArray(
+        sc.array(
+            dims=['x'],
+            values=_quadratic(x.values),
+            variances=2.0 + x.values,
+            unit='counts',
+        ),
+        coords={'x': x},
+    )
+
+    with pytest.raises(sc.VariancesError, match="signals with variances"):
+        smooth_relative(data, scale=0.1)
+
+
+@pytest.mark.parametrize(
+    ("x", "y"),
+    [
+        (np.arange(1.0, 4.0), sc.ones(dims=['x'], shape=[3])),
+        (sc.arange('x', 1.0, 4.0), np.ones(3)),
+    ],
+)
+def test_rejects_non_scipp_pair_member(x, y):
+    with pytest.raises(TypeError, match="DataArray or a pair of Variables"):
+        smooth_relative(x, y, scale=0.1)
+
+
+def test_rejects_variables_with_different_dimensions():
+    x = sc.arange('x', 1.0, 4.0)
+    y = sc.ones(dims=['y'], shape=[3])
+
+    with pytest.raises(sc.DimensionError, match="same dimension"):
+        smooth_relative(x, y, scale=0.1)
+
+
+def test_rejects_data_array_without_dimension_coordinate():
+    data = sc.DataArray(sc.ones(dims=['x'], shape=[3]))
+
+    with pytest.raises(sc.CoordError, match="dimension coordinate"):
+        smooth_relative(data, scale=0.1)
+
+
+def test_rejects_data_array_with_bin_edge_coordinate():
+    data = sc.DataArray(
+        sc.ones(dims=['x'], shape=[3]),
+        coords={'x': sc.arange('x', 1.0, 5.0)},
+    )
+
+    with pytest.raises(sc.CoordError, match="bin edges"):
+        smooth_relative(data, scale=0.1)
+
+
+def test_rejects_data_array_with_masks():
+    data = sc.DataArray(
+        sc.ones(dims=['x'], shape=[3]),
+        coords={'x': sc.arange('x', 1.0, 4.0)},
+        masks={'bad': sc.array(dims=['x'], values=[False, True, False])},
+    )
+
+    with pytest.raises(ValueError, match="data with masks"):
+        smooth_relative(data, scale=0.1)
+
+
+def test_rejects_y_with_data_array():
+    x, y = _variables([1.0, 2.0], [3.0, 4.0])
+    data = sc.DataArray(y, coords={'x': x})
+
+    with pytest.raises(TypeError, match="y must be omitted"):
+        smooth_relative(data, y, scale=0.1)
+
+
+def test_rejects_geometric_grid_larger_than_limit():
+    x = np.array([1.0, 1.01, 2.0])
+    y = np.ones_like(x)
+
+    with pytest.raises(
+        ValueError,
+        match=r"geometric resampling would require too many points.*max_grid_points=70",
+    ):
+        _smooth_values(smooth_relative, x, y, scale=0.1, max_grid_points=70)
+
+
+def test_accepts_geometric_grid_equal_to_limit():
+    x = np.array([1.0, 1.01, 2.0])
+    y = np.ones_like(x)
+
+    actual = _smooth_values(smooth_relative, x, y, scale=0.1, max_grid_points=71)
+
+    np.testing.assert_allclose(actual, y)
+
+
+def test_geometric_input_does_not_gain_a_point_from_roundoff():
+    size = 100
+    x = _geometric_cell_centers(0.1, 0.9, size)
+    y = _quadratic(x)
+
+    actual = _smooth_values(
+        smooth_relative,
+        x,
+        y,
+        scale=0.1,
+        max_grid_points=size,
+    )
+
+    assert actual.shape == y.shape
+
+
+@pytest.mark.filterwarnings("ignore:divide by zero encountered in scalar divide")
+@pytest.mark.parametrize("x0", [1.0, 3.0, 10.0, 0.1])
+def test_pathologically_close_coordinates_fail_before_allocation(x0):
+    x = np.array([x0, np.nextafter(x0, 2.0 * x0), 2.0 * x0])
+    y = np.ones_like(x)
+
+    with pytest.raises(ValueError, match="exceeding max_grid_points=1,000,000"):
+        _smooth_values(smooth_relative, x, y, scale=0.1)
+
+
+def test_wide_coordinate_range_does_not_overflow_grid_construction():
+    x = np.array([1e-300, 1.0, 1e300])
+    y = np.ones_like(x)
+
+    actual = _smooth_values(smooth_relative, x, y, scale=0.1)
+
+    np.testing.assert_allclose(actual, y)
+
+
+def test_kernel_stencil_is_bounded_before_allocation():
+    # A gaussian truncated at tail=1e-12 spans roughly 14 sigma, which is about
+    # 1e7 offsets on this grid. The stencil must be clamped to offsets that can
+    # reach the input rather than allocated at its nominal width.
+    size = 21
+    x = np.geomspace(1.0, np.exp((size - 1) * 1e-6), size)
+    y = np.arange(size, dtype=float)
+
+    actual = _smooth_values(smooth_relative, x, y, scale=1.0)
+
+    # The kernel is flat to within 1e-5 over the whole input, so every point
+    # averages the entire array.
+    np.testing.assert_allclose(actual, np.full(size, y.mean()), rtol=1e-4)
+
+
+def test_asymmetric_kernel_convolution_matches_direct_weighted_sum():
+    size = 20
+    log_spacing = 0.03
+    alpha = 0.2
+    kernel = uniform(loc=0.5, scale=1.0)
+    y = np.arange(size, dtype=float) ** 2
+    x = np.geomspace(1.0, np.exp((size - 1) * log_spacing), size)
+
+    offsets = np.arange(-size + 1, size)
+    lower = np.expm1((offsets - 0.5) * log_spacing) / alpha
+    upper = np.expm1((offsets + 0.5) * log_spacing) / alpha
+    weights = kernel.cdf(upper) - kernel.cdf(lower)
+    expected = np.empty_like(y)
+    for i in range(size):
+        valid = (0 <= i + offsets) & (i + offsets < size)
+        denominator = np.sum(weights[valid])
+        expected[i] = (
+            np.dot(weights[valid], y[i + offsets[valid]]) / denominator
+            if denominator > 0.0
+            else np.nan
+        )
+
+    actual = _smooth_values(
+        smooth_relative,
+        x,
+        y,
+        scale=alpha,
+        kernel=kernel,
+        tail=1e-12,
+    )
+
+    np.testing.assert_allclose(actual, expected, equal_nan=True)
+
+
+def test_kernel_with_no_reachable_mass_returns_nan():
+    y = np.arange(5.0)
+    x = np.geomspace(1.0, np.exp((y.size - 1) * 0.1), y.size)
+    actual = _smooth_values(
+        smooth_relative,
+        x,
+        y,
+        scale=0.1,
+        kernel=uniform(loc=100.0, scale=1.0),
+        tail=1e-12,
+    )
+
+    assert np.all(np.isnan(actual))
+
+
+@pytest.mark.parametrize("scale", [np.nan, np.inf, -np.inf, -1.0])
+def test_rejects_invalid_relative_scale_before_noop_return(scale):
+    x, y = _variables([], [])
+    with pytest.raises(ValueError, match="scale must be non-negative"):
+        smooth_relative(x, y, scale=scale)
+
+
+@pytest.mark.parametrize("tail", [np.nan, np.inf, -np.inf, 0.0, 1.0])
+def test_rejects_invalid_tail_before_noop_return(tail):
+    x, y = _variables([], [])
+    with pytest.raises(ValueError, match="tail must be between 0 and 1"):
+        smooth_relative(x, y, scale=0.1, tail=tail)
+
+
+@pytest.mark.parametrize(
+    ("x", "message"),
+    [
+        ([-1.0], "x must be positive"),
+        ([np.nan], "x must contain only finite values"),
+        ([np.inf], "x must contain only finite values"),
+    ],
+)
+def test_rejects_invalid_single_coordinate_before_noop_return(x, message):
+    x, y = _variables(x, [1.0])
+    with pytest.raises(ValueError, match=message):
+        smooth_relative(x, y, scale=0.1)
+
+
+@pytest.mark.parametrize(
+    ("function", "kwargs"),
+    [
+        (smooth_relative, {"scale": 0.1}),
+        (smooth, {"scale": sc.scalar(0.1, unit='m')}),
+    ],
+)
+@pytest.mark.parametrize("max_grid_points", [True, 2.5])
+def test_rejects_non_integer_max_grid_points(function, kwargs, max_grid_points):
+    x, y = _variables([], [])
+    with pytest.raises(TypeError, match="max_grid_points must be an integer"):
+        function(x, y, max_grid_points=max_grid_points, **kwargs)
+
+
+@pytest.mark.parametrize(
+    ("function", "kwargs"),
+    [
+        (smooth_relative, {"scale": 0.1}),
+        (smooth, {"scale": sc.scalar(0.1, unit='m')}),
+    ],
+)
+@pytest.mark.parametrize("max_grid_points", [-1, 0, 1])
+def test_rejects_too_small_max_grid_points(function, kwargs, max_grid_points):
+    x, y = _variables([], [])
+    with pytest.raises(ValueError, match="max_grid_points must be at least 2"):
+        function(x, y, max_grid_points=max_grid_points, **kwargs)