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..6dedaf603
--- /dev/null
+++ b/src/scippneutron/smoothing.py
@@ -0,0 +1,620 @@
+# 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
+
+import warnings
+from numbers import Real
+from typing import Any, Protocol, 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]: ...
+
+
+type _Kernel = str | _Distribution
+type _FloatArray = NDArray[np.float64]
+type _IntArray = NDArray[np.int64]
+
+_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[:0], weights[:0]
+
+    first = nonzero[0]
+    last = nonzero[-1] + 1
+    return offsets[first:last], weights[first:last]
+
+
+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)
+    support_min, support_max = _kernel_support(dist)
+
+    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")
+    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))
+        )
+        # The geometric domain boundary may leave u_left non-finite; m_min is
+        # safely clamped to the finite input below. u_right has no such bound.
+        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.
+
+    Sum weights whose offsets remain within the input at each output position.
+    This uses cumulative sums and is 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:])
+
+    out = np.zeros(n, dtype=float)
+    valid = lower <= upper
+    out[valid] = cumsum[upper[valid] - m_min + 1] - cumsum[lower[valid] - m_min]
+    return out
+
+
+def _smooth_with_weights(
+    y: _FloatArray, offsets: _IntArray, weights: _FloatArray
+) -> _FloatArray:
+    if weights.size == 0:
+        return np.full_like(y, np.nan)
+
+    # 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])
+    source_begin = max(start, 0)
+    source_end = min(start + y.size, full.size)
+    destination_begin = source_begin - start
+    destination_end = destination_begin + source_end - source_begin
+    numerator = np.zeros_like(y)
+    numerator[destination_begin:destination_end] = full[source_begin:source_end]
+    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 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 _validate_data(data: object) -> sc.Variable:
+    """Validate user input and return its dimension coordinate."""
+    if not isinstance(data, sc.DataArray):
+        raise TypeError("expected a DataArray")
+    if data.ndim != 1:
+        raise sc.DimensionError("data must be one-dimensional")
+    if data.is_binned:
+        raise sc.DTypeError("data must not be binned")
+    if data.dim not in data.coords:
+        raise sc.CoordError("data must have a dimension coordinate")
+    coord = data.coords[data.dim]
+    if coord.sizes != data.sizes:
+        raise sc.CoordError(
+            "the dimension coordinate must have the same shape as the data"
+        )
+    if data.masks:
+        raise ValueError("smoothing data with masks is not supported")
+    if data.variances is not None:
+        raise sc.VariancesError(
+            "Smoothing signals with variances is not supported because it would "
+            "introduce correlations between data points."
+        )
+    if np.any(~np.isfinite(data.values)):
+        warnings.warn(
+            "Data contains NaNs or infinities; smoothing may fall back to a slower "
+            "method.",
+            UserWarning,
+            stacklevel=3,
+        )
+    return coord
+
+
+def _with_values(data: sc.DataArray, values: _FloatArray) -> sc.DataArray:
+    # The new container shares unchanged coordinates, but its data is independent.
+    out = data.copy(deep=False)
+    out.data = sc.array(dims=data.dims, values=values, unit=data.unit)
+    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(
+    data: sc.DataArray,
+    *,
+    scale: sc.Variable,
+    kernel: _Kernel = "gaussian",
+    tail: float = 1e-12,
+    max_grid_points: int = 1_000_000,
+) -> sc.DataArray:
+    """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
+    ----------
+    data:
+        One-dimensional data to smooth. Must have a strictly increasing
+        dimension coordinate.
+    scale:
+        Finite, non-negative 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. Must be strictly between zero and one.
+    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. Coordinates and units are preserved. Points where no
+        kernel mass falls inside the input domain are NaN; this can occur at
+        the boundaries for one-sided kernels.
+
+    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 input is not one-dimensional or ``scale`` is not scalar.
+    scipp.CoordError
+        If a data array has no dimension coordinate or has a bin-edge
+        coordinate.
+    scipp.DTypeError
+        If ``data`` is binned.
+    scipp.UnitError
+        If the unit of ``scale`` is incompatible with the coordinate unit.
+    scipp.VariancesError
+        If the signal or ``scale`` has variances.
+    TypeError
+        If ``data`` is not a data array, ``scale`` is not a variable, ``kernel``
+        is not a distribution-like object, or ``max_grid_points`` is not an
+        integer.
+
+    Warns
+    -----
+    UserWarning
+        If the data contains NaNs or infinities, since smoothing may fall back to
+        a slower method.
+    """
+    x = _validate_data(data)
+    values = _smooth_values(
+        _Uniform(),
+        x.values,
+        data.values,
+        scale=_scale_in_coordinate_unit(scale, x),
+        kernel=kernel,
+        tail=tail,
+        max_grid_points=max_grid_points,
+    )
+    return _with_values(data, values)
+
+
+def smooth_relative(
+    data: sc.DataArray,
+    *,
+    scale: float | sc.Variable,
+    kernel: _Kernel = "gaussian",
+    tail: float = 1e-12,
+    max_grid_points: int = 1_000_000,
+) -> sc.DataArray:
+    """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
+    ----------
+    data:
+        One-dimensional data to smooth. Must have a positive, strictly
+        increasing dimension coordinate.
+    scale:
+        Finite, non-negative, 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. Must be strictly
+        between zero and one.
+    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. Coordinates and units are preserved. Points where no
+        kernel mass falls inside the input domain are NaN; this can occur at
+        the boundaries for one-sided kernels.
+
+    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 input is not one-dimensional or ``scale`` is not scalar.
+    scipp.CoordError
+        If a data array has no dimension coordinate or has a bin-edge
+        coordinate.
+    scipp.DTypeError
+        If ``data`` is binned.
+    scipp.UnitError
+        If ``scale`` is a variable with a non-dimensionless unit.
+    scipp.VariancesError
+        If the signal or ``scale`` has variances.
+    TypeError
+        If ``data`` is not a data array, ``scale`` is neither a real number nor a
+        variable, ``kernel`` is not a distribution-like object, or
+        ``max_grid_points`` is not an integer.
+
+    Warns
+    -----
+    UserWarning
+        If the data contains NaNs or infinities, since smoothing may fall back to
+        a slower method.
+    """
+    x = _validate_data(data)
+    values = _smooth_values(
+        _Geometric(),
+        x.values,
+        data.values,
+        scale=_dimensionless_scale(scale),
+        kernel=kernel,
+        tail=tail,
+        max_grid_points=max_grid_points,
+    )
+    return _with_values(data, values)
diff --git a/tests/smoothing_test.py b/tests/smoothing_test.py
new file mode 100644
index 000000000..e5cfe6af3
--- /dev/null
+++ b/tests/smoothing_test.py
@@ -0,0 +1,904 @@
+# 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 _data_array(x, y):
+    return sc.DataArray(
+        sc.array(
+            dims=['x'],
+            values=np.asarray(y),
+            unit='counts',
+        ),
+        coords={'x': sc.array(dims=['x'], values=np.asarray(x), unit='m')},
+    )
+
+
+def _quadratic_smoothing_error(*, size, alpha, tail=1e-9):
+    lower = 0.1
+    upper = 0.9
+    x = _geometric_cell_centers(lower, upper, size)
+    actual = smooth_relative(
+        _data_array(x, _quadratic(x)),
+        scale=alpha,
+        tail=tail,
+    ).values
+    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)
+    actual = smooth(
+        _data_array(x, _quadratic(x)),
+        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)
+    actual = smooth(
+        _data_array(x, _quadratic(x)),
+        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')
+    data = sc.DataArray(y, coords={'x': x})
+
+    in_meters = smooth(data, scale=sc.scalar(0.1, unit='m'))
+    in_centimeters = smooth(data, 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')
+    data = sc.DataArray(y, coords={'x': x})
+    scale = sc.scalar(0.1, unit='m')
+
+    actual = smooth(
+        data,
+        scale=scale,
+        kernel=uniform(loc=-1.0, scale=2.0),
+    )
+    expected = smooth(data, 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')
+    data = sc.DataArray(y, coords={'x': x})
+
+    actual = smooth(
+        data,
+        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
+    data = sc.DataArray(sc.array(dims=['x'], values=values), coords={'x': x})
+
+    with pytest.warns(UserWarning, match="NaNs or infinities"):
+        actual = smooth(data, 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_zero_weight_does_not_propagate_nonfinite_value():
+    data = _data_array([0.0, 1.0, 2.0], [np.inf, 2.0, 3.0])
+
+    with pytest.warns(UserWarning, match="NaNs or infinities"):
+        actual = smooth(
+            data,
+            scale=sc.scalar(1.0, unit='m'),
+            kernel=uniform(loc=0.5, scale=1.0),
+        )
+
+    np.testing.assert_equal(actual.values, [2.0, 3.0, np.nan])
+
+
+@pytest.mark.parametrize(
+    ("function", "x", "scale"),
+    [
+        (smooth, [0.0, 1.0], sc.scalar(0.1, unit='m')),
+        (smooth_relative, [1.0, 2.0], 0.1),
+    ],
+)
+@pytest.mark.parametrize("nonfinite", [np.nan, np.inf, -np.inf])
+def test_warns_for_nonfinite_data(function, x, scale, nonfinite):
+    data = _data_array(x, [1.0, nonfinite])
+
+    with pytest.warns(UserWarning, match="NaNs or infinities.*slower") as emitted:
+        function(data, scale=scale)
+
+    assert emitted[0].filename == __file__
+
+
+def test_fixed_width_smoothing_preserves_coordinates():
+    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'))
+    assert actual.sizes == data.sizes
+    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')
+
+    actual = smooth(_data_array(x_values, y_values), scale=scale)
+
+    xp_values = np.linspace(0.0, 1.0, 5)
+    smoothed_yp = smooth(
+        _data_array(xp_values, np.interp(xp_values, x_values, y_values)), 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))
+    actual = smooth(
+        _data_array(x, _quadratic(x)),
+        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():
+    data = _data_array([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(data, scale=sc.scalar(0.1, unit='m'), max_grid_points=4)
+
+
+def test_fixed_width_smoothing_accepts_resampled_grid_equal_to_limit():
+    data = _data_array([0.0, 0.25, 1.0], [1.0, 1.0, 1.0])
+
+    actual = smooth(
+        data,
+        scale=sc.scalar(0.1, unit='m'),
+        max_grid_points=5,
+    )
+
+    np.testing.assert_allclose(actual.values, data.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 = function(_data_array(x, y), max_grid_points=10, **kwargs).values
+
+    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 = function(_data_array(x, y), max_grid_points=4, **kwargs).values
+
+    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,
+):
+    data = _data_array(coordinates, [1.0, 1.0, 1.0])
+
+    with pytest.raises(ValueError, match="exceeding max_grid_points=1,000,000"):
+        smooth(data, 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):
+    data = _data_array([], [])
+
+    with pytest.raises(ValueError, match="scale must be non-negative"):
+        smooth(data, scale=sc.scalar(scale, unit='m'))
+
+
+def test_fixed_width_smoothing_rejects_non_scalar_scale():
+    data = _data_array([0.0, 1.0], [1.0, 2.0])
+
+    with pytest.raises(sc.DimensionError, match="scale must be a scalar"):
+        smooth(data, scale=sc.array(dims=['scale'], values=[0.1], unit='m'))
+
+
+def test_fixed_width_smoothing_rejects_non_variable_scale():
+    data = _data_array([0.0, 1.0], [1.0, 2.0])
+
+    with pytest.raises(TypeError, match=r"scale must be a scipp\.Variable"):
+        smooth(data, scale=0.1)  # type: ignore[arg-type]
+
+
+def test_fixed_width_smoothing_rejects_incompatible_scale_unit():
+    data = _data_array([0.0, 1.0], [1.0, 2.0])
+
+    with pytest.raises(sc.UnitError):
+        smooth(data, scale=sc.scalar(0.1, unit='s'))
+
+
+def test_fixed_width_smoothing_rejects_scale_with_variance():
+    data = _data_array([0.0, 1.0], [1.0, 2.0])
+
+    with pytest.raises(sc.VariancesError, match="scales with variances"):
+        smooth(
+            data,
+            scale=sc.scalar(0.1, variance=0.01, unit='m'),
+        )
+
+
+def test_relative_width_gaussian_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_relative_width_gaussian_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_relative_width_gaussian_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_relative_width_gaussian_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_relative_width_compact_kernel_matches_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_relative(
+        _data_array(x, _quadratic(x)),
+        scale=alpha,
+        kernel=kernel,
+    ).values
+    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
+
+
+@pytest.mark.parametrize(
+    ("function", "x", "scale"),
+    [
+        (smooth, [0.0, 1.0], sc.scalar(0.0, unit='m')),
+        (smooth_relative, [1.0, 2.0], 0.0),
+    ],
+)
+def test_zero_scale_returns_independent_copy(function, x, scale):
+    data = _data_array(x, [1.0, 2.0])
+
+    actual = function(data, scale=scale)
+
+    assert sc.identical(actual, data)
+    actual.values[0] = -1.0
+    assert data.values[0] == 1.0
+
+
+@pytest.mark.parametrize(
+    ("function", "scale"),
+    [
+        (smooth, sc.scalar(0.1, unit='m')),
+        (smooth_relative, 0.1),
+    ],
+)
+def test_rejects_unknown_kernel_name(function, scale):
+    data = _data_array([1.0, 2.0], [1.0, 2.0])
+
+    with pytest.raises(ValueError, match="unknown kernel 'unknown'"):
+        function(data, scale=scale, kernel="unknown")
+
+
+@pytest.mark.parametrize(
+    ("function", "scale"),
+    [
+        (smooth, sc.scalar(0.1, unit='m')),
+        (smooth_relative, 0.1),
+    ],
+)
+def test_rejects_non_monotonic_coordinate(function, scale):
+    data = _data_array([1.0, 3.0, 2.0], [1.0, 2.0, 3.0])
+
+    with pytest.raises(ValueError, match="x must be strictly increasing"):
+        function(data, scale=scale)
+
+
+@pytest.mark.parametrize(
+    ("function", "scale"),
+    [
+        (smooth, sc.scalar(0.1, unit='m')),
+        (smooth_relative, 0.1),
+    ],
+)
+def test_rejects_two_dimensional_input(function, scale):
+    data = sc.DataArray(sc.ones(dims=['x', 'y'], shape=[2, 2]))
+
+    with pytest.raises(sc.DimensionError, match="data must be one-dimensional"):
+        function(data, scale=scale)
+
+
+def test_relative_smoothing_preserves_metadata_and_input():
+    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)
+
+    assert sc.identical(data, original)
+    assert isinstance(actual, sc.DataArray)
+    assert actual.sizes == data.sizes
+    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():
+    data = _data_array([1.0, 2.0, 3.0], [1.0, 4.0, 9.0])
+
+    actual = smooth_relative(data, scale=sc.scalar(0.1))
+    expected = smooth_relative(data, scale=0.1)
+
+    assert sc.identical(actual, expected)
+
+
+def test_relative_smoothing_rejects_non_scalar_scale():
+    data = _data_array([1.0, 2.0], [1.0, 2.0])
+
+    with pytest.raises(sc.DimensionError, match="scale must be a scalar"):
+        smooth_relative(data, scale=sc.array(dims=['scale'], values=[0.1]))
+
+
+def test_relative_smoothing_rejects_scale_with_unit():
+    data = _data_array([1.0, 2.0], [1.0, 2.0])
+
+    with pytest.raises(sc.UnitError):
+        smooth_relative(data, scale=sc.scalar(0.1, unit='m'))
+
+
+def test_relative_smoothing_rejects_scale_with_variance():
+    data = _data_array([1.0, 2.0], [1.0, 2.0])
+
+    with pytest.raises(sc.VariancesError, match="scales with variances"):
+        smooth_relative(data, scale=sc.scalar(0.1, variance=0.01))
+
+
+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(
+    ("function", "scale"),
+    [
+        (smooth, sc.scalar(0.1, unit='m')),
+        (smooth_relative, 0.1),
+    ],
+)
+def test_rejects_non_data_array(function, scale):
+    with pytest.raises(TypeError, match="expected a DataArray"):
+        function(sc.arange('x', 1.0, 4.0), scale=scale)
+
+
+@pytest.mark.parametrize(
+    ("function", "scale"),
+    [
+        (smooth, sc.scalar(0.1, unit='m')),
+        (smooth_relative, 0.1),
+    ],
+)
+def test_rejects_binned_data(function, scale):
+    binned = sc.bins(
+        dim='event',
+        data=sc.arange('event', 4.0),
+        begin=sc.array(dims=['x'], values=[0, 2], unit=None),
+        end=sc.array(dims=['x'], values=[2, 4], unit=None),
+    )
+    data = sc.DataArray(
+        binned, coords={'x': sc.array(dims=['x'], values=[1.0, 2.0], unit='m')}
+    )
+
+    with pytest.raises(sc.DTypeError, match="data must not be binned"):
+        function(data, scale=scale)
+
+
+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_scalar_dimension_coordinate():
+    data = sc.DataArray(sc.ones(dims=['x'], shape=[3]), coords={'x': sc.scalar(1.0)})
+
+    with pytest.raises(sc.CoordError, match="same shape as the data"):
+        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="same shape as the data"):
+        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_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_relative(_data_array(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_relative(_data_array(x, y), scale=0.1, max_grid_points=71).values
+
+    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_relative(
+        _data_array(x, y),
+        scale=0.1,
+        max_grid_points=size,
+    ).values
+
+    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_relative(_data_array(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_relative(_data_array(x, y), scale=0.1).values
+
+    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_relative(_data_array(x, y), scale=1.0).values
+
+    # 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)
+
+
+@pytest.mark.parametrize(
+    "kernel",
+    [
+        pytest.param(uniform(loc=0.5, scale=1.0), id="positive-offsets"),
+        pytest.param(uniform(loc=-1.5, scale=1.0), id="negative-offsets"),
+    ],
+)
+def test_asymmetric_kernel_convolution_matches_direct_weighted_sum(kernel):
+    size = 20
+    log_spacing = 0.03
+    alpha = 0.2
+    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_relative(
+        _data_array(x, y),
+        scale=alpha,
+        kernel=kernel,
+        tail=1e-12,
+    ).values
+
+    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_relative(
+        _data_array(x, y),
+        scale=0.1,
+        kernel=uniform(loc=100.0, scale=1.0),
+        tail=1e-12,
+    ).values
+
+    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):
+    data = _data_array([], [])
+    with pytest.raises(ValueError, match="scale must be non-negative"):
+        smooth_relative(data, 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):
+    data = _data_array([], [])
+    with pytest.raises(ValueError, match="tail must be between 0 and 1"):
+        smooth_relative(data, 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):
+    data = _data_array(x, [1.0])
+    with pytest.raises(ValueError, match=message):
+        smooth_relative(data, 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):
+    data = _data_array([], [])
+    with pytest.raises(TypeError, match="max_grid_points must be an integer"):
+        function(data, 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):
+    data = _data_array([], [])
+    with pytest.raises(ValueError, match="max_grid_points must be at least 2"):
+        function(data, max_grid_points=max_grid_points, **kwargs)