diff --git a/README.md b/README.md
index bc474098..2b8b49d8 100644
--- a/README.md
+++ b/README.md
@@ -125,7 +125,7 @@ result, and prediction interfaces.
 
 | Formula surface    | Supported terms                                                                                      |
 | ------------------ | ---------------------------------------------------------------------------------------------------- |
-| Univariate smooths | `s(..., bs='cr')`, `cs`, `cc`, `cp`, `ps`, `tp`, `ts`                                                |
+| Univariate smooths | `s(..., bs='bs')`, `cr`, `cs`, `cc`, `cp`, `ps`, `tp`, `ts`                                          |
 | Structured smooths | random effects `re`, factor smooths `fs`, sum-to-zero factor smooths `sz`                            |
 | Tensor products    | `te(...)` and `ti(...)` over supported numeric marginals                                             |
 | Parametric terms   | numeric and factor terms, supported interactions, intercept policies, and formula offsets            |
diff --git a/nampy/gam/compiler/factory.py b/nampy/gam/compiler/factory.py
index f82708fa..c40cb0fc 100644
--- a/nampy/gam/compiler/factory.py
+++ b/nampy/gam/compiler/factory.py
@@ -15,6 +15,7 @@
 from ..smooths.registry import make_smooth_term
 from ..smooths.shape.bivariate import BivariateShapePSplineTerm
 from ..smooths.shape.scop import ShapeConstrainedPSplineTerm
+from ..smooths.univariate.bs import DerivativeBSplineTerm1D
 from ..smooths.univariate.cr import CubicSplineTerm
 from ..smooths.univariate.ps import PSplineTerm1D
 from ..specs import LinearPredictorSpec, PenaltyGroupSpec, TermSpec
@@ -22,6 +23,7 @@
     CubicRegressionSmoothSpec,
     CubicShrinkageSmoothSpec,
     CyclicCubicRegressionSmoothSpec,
+    DerivativeBSplineSmoothSpec,
     FactorSmoothInteractionSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
@@ -162,6 +164,28 @@ def instantiate_term(term_like: TermSpec | Any):
             metadata=metadata,
         )
 
+    if isinstance(smooth_spec, DerivativeBSplineSmoothSpec):
+        if len(features) != 1:
+            raise NotImplementedError(
+                "Current runtime only materializes 1D s(..., bs='bs') terms."
+            )
+        return DerivativeBSplineTerm1D(
+            feature=features[0],
+            k=smooth_spec.k,
+            m=smooth_spec.m,
+            label=label,
+            term_id=term_like.term_id,
+            smoothing_id=smoothing_id,
+            by=by,
+            sp=smooth_spec.sp,
+            select=smooth_spec.select,
+            fixed=smooth_spec.fx,
+            constraint_mode=smooth_spec.constraint_mode,
+            pc=smooth_spec.pc,
+            knots=smooth_spec.knots,
+            metadata=metadata,
+        )
+
     if isinstance(smooth_spec, ShapeConstrainedSmoothSpec):
         if len(features) == 2:
             return BivariateShapePSplineTerm(
@@ -242,6 +266,7 @@ def instantiate_term(term_like: TermSpec | Any):
             by=by,
             sp=smooth_spec.sp,
             select=smooth_spec.select,
+            m=smooth_spec.m,
             xt=smooth_spec.xt,
             fixed=smooth_spec.fx,
             knots=smooth_spec.knots,
@@ -258,6 +283,7 @@ def instantiate_term(term_like: TermSpec | Any):
             by=by,
             sp=smooth_spec.sp,
             select=smooth_spec.select,
+            m=smooth_spec.m,
             xt=smooth_spec.xt,
             fixed=smooth_spec.fx,
             knots=smooth_spec.knots,
@@ -304,9 +330,15 @@ def _expected_penalty_group_size(runtime_term):
     if bool(getattr(runtime_term, "fixed", False)):
         return 0
 
+    if hasattr(runtime_term, "expected_linked_penalty_count"):
+        value = runtime_term.expected_linked_penalty_count
+        return None if value is None else int(value)
+
     fixed_flags = getattr(runtime_term, "fixed_flags", None)
     if fixed_flags is not None:
         n_penalties = int(np.sum(~np.asarray(fixed_flags, dtype=bool)))
+    elif getattr(runtime_term, "n_main_penalties", None) is not None:
+        n_penalties = int(runtime_term.n_main_penalties)
     else:
         term_type = str(getattr(runtime_term, "term_type", "smooth"))
         if term_type in {"tensor_smooth", "tensor_interaction"}:
diff --git a/nampy/gam/smooths/__init__.py b/nampy/gam/smooths/__init__.py
index 1fc77f41..5c86567c 100644
--- a/nampy/gam/smooths/__init__.py
+++ b/nampy/gam/smooths/__init__.py
@@ -19,6 +19,7 @@
 )
 from .tensor.te import TensorProductSplineTerm
 from .tensor.ti import InteractionTensorProductSplineTerm
+from .univariate.bs import DerivativeBSplineTerm1D
 from .univariate.cr import CubicSplineTerm
 from .univariate.ps import PSplineTerm1D
 from .univariate.tp import ThinPlateSplineTerm
@@ -28,6 +29,7 @@
 te = TensorProductSplineTerm
 ti = InteractionTensorProductSplineTerm
 
+bs = DerivativeBSplineTerm1D
 cr = cs = cc = CubicSplineTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
@@ -55,6 +57,7 @@
     "sync_by_state_attributes",
     "build_penalty_definition",
     "CubicSplineTerm",
+    "DerivativeBSplineTerm1D",
     "PSplineTerm1D",
     "ThinPlateSplineTerm",
     "TensorProductSplineTerm",
@@ -65,6 +68,7 @@
     "ShapeConstrainedPSplineTerm",
     "te",
     "ti",
+    "bs",
     "cr",
     "cs",
     "cc",
diff --git a/nampy/gam/smooths/categorical/fs.py b/nampy/gam/smooths/categorical/fs.py
index 414e455f..3bcf2e4b 100644
--- a/nampy/gam/smooths/categorical/fs.py
+++ b/nampy/gam/smooths/categorical/fs.py
@@ -15,6 +15,7 @@
 from ..algebra import rowwise_kronecker
 from ..registry import make_smooth_term
 from ..smooth_base import BaseSmoothTerm, by_values_from_new_data, column_as_object
+from ..univariate.bs import DerivativeBSplineTerm1D
 from ..univariate.cr import CubicSplineTerm
 from ..univariate.ps import PSplineTerm1D
 from .categorical_utils import (
@@ -97,6 +98,7 @@ def _build_base_smooth_term(
     by,
     knots,
     xt_rest,
+    outer_m,
     mode,  # "fs" or "sz"
     select,
     constraint_mode,
@@ -106,7 +108,7 @@ def _build_base_smooth_term(
     Build the per-level base smooth used inside fs/sz.
 
     Supported base smooth classes in the current codebase:
-    cr, cs, cc, cp, ps, tp, ts
+    bs, cr, cs, cc, cp, ps, tp, ts
     """
     base_bs = str(base_bs).lower()
     metric_features = list(metric_features)
@@ -123,9 +125,9 @@ def _build_base_smooth_term(
             f"for bs in {{'tp','ts'}}, got base bs={base_bs!r}."
         )
 
-    if xt_rest is not None and base_bs not in {"tp", "ts", "ps", "cp"}:
+    if xt_rest is not None and base_bs not in {"bs", "tp", "ts", "ps", "cp"}:
         raise NotImplementedError(
-            "Extra xt options are currently only supported for tp/ts/ps/cp base "
+            "Extra xt options are currently only supported for bs/tp/ts/ps/cp base "
             "smooths, "
             f"got xt={xt_rest!r} with base bs={base_bs!r}."
         )
@@ -149,7 +151,9 @@ def _build_base_smooth_term(
         )
 
     if base_bs in {"ps", "cp"}:
-        ps_m = None if xt_rest is None else xt_rest.get("m", None)
+        ps_m = outer_m
+        if ps_m is None and xt_rest is not None:
+            ps_m = xt_rest.get("m", None)
         # For fs/sz, mgcv keeps the outer basis dimension and uses xt mainly to
         # choose the base smoother family / order parameters.
         ps_k = k
@@ -170,6 +174,26 @@ def _build_base_smooth_term(
             metadata=metadata,
         )
 
+    if base_bs == "bs":
+        bs_m = outer_m
+        if bs_m is None and xt_rest is not None:
+            bs_m = xt_rest.get("m", None)
+        return DerivativeBSplineTerm1D(
+            feature=metric_features[0],
+            k=k,
+            m=bs_m,
+            label=label,
+            smoothing_id=None,
+            by=by,
+            sp=None,
+            select=bool(select),
+            fixed=bool(fixed),
+            constraint_mode=str(constraint_mode),
+            pc=None,
+            knots=knots,
+            metadata=metadata,
+        )
+
     if base_bs in {"tp", "ts"}:
         return make_smooth_term(
             base_bs,
@@ -192,11 +216,13 @@ def _build_base_smooth_term(
 
     raise NotImplementedError(
         f"Current {mode} implementation supports base bs in "
-        f"{{'cr','cs','cc','cp','ps','tp','ts'}}, got {base_bs!r}."
+        f"{{'bs','cr','cs','cc','cp','ps','tp','ts'}}, got {base_bs!r}."
     )
 
 
 def _penalty_rank_from_base_term(base_term, basis_matrix, penalty_matrix) -> int:
+    if isinstance(base_term, DerivativeBSplineTerm1D):
+        return int(base_term._setup.ranks[0])
     if isinstance(base_term, PSplineTerm1D) and len(base_term.penalties) > 0:
         if str(base_term.basis_name).lower() == "cp":
             return int(base_term._setup.rank)
@@ -286,6 +312,7 @@ def __init__(
         by=None,
         sp=None,
         select=False,
+        m=None,
         xt=None,
         fixed=False,
         knots=None,
@@ -307,6 +334,7 @@ def __init__(
         self.term_type = term_type
         self.k = int(k)
         self.select = bool(select)
+        self.m = m
         self.xt = xt
         self.fixed = bool(fixed)
         self.knots = knots
@@ -438,6 +466,7 @@ def _build_delegate_base_or_re(self, X, feature_names, *, default_bs, mode):
                 by=self.by,
                 knots=self.knots,
                 xt_rest=base_spec.xt_rest,
+                outer_m=self.m,
                 mode=mode,
                 select=self.select,
                 constraint_mode=("auto" if mode == "fs" else "never"),
@@ -504,6 +533,7 @@ def __init__(
         by=None,
         sp=None,
         select=False,
+        m=None,
         xt=None,
         fixed=False,
         knots=None,
@@ -520,6 +550,7 @@ def __init__(
             by=by,
             sp=sp,
             select=select,
+            m=m,
             xt=xt,
             fixed=fixed,
             knots=knots,
@@ -558,6 +589,7 @@ def fit(self, X, feature_names):
             by=None,
             knots=self.knots,
             xt_rest=base_spec.xt_rest,
+            outer_m=self.m,
             mode="fs",
             select=False,
             constraint_mode=base_constraint,
@@ -567,7 +599,7 @@ def fit(self, X, feature_names):
 
         if len(base_term.penalties) > 1:
             raise NotImplementedError(
-                'bs="fs" currently requires a singly penalized base smooth.'
+                '"fs" smooth cannot use a multiply penalized basis (wrong basis in xt)'
             )
 
         self._base_term = base_term
@@ -789,6 +821,7 @@ def __init__(
         by=None,
         sp=None,
         select=False,
+        m=None,
         xt=None,
         fixed=False,
         knots=None,
@@ -805,6 +838,7 @@ def __init__(
             by=by,
             sp=sp,
             select=select,
+            m=m,
             xt=xt,
             fixed=fixed,
             knots=knots,
@@ -834,6 +868,7 @@ def fit(self, X, feature_names):
             by=None,
             knots=self.knots,
             xt_rest=base_spec.xt_rest,
+            outer_m=self.m,
             mode="sz",
             select=False,
             constraint_mode="never",
@@ -843,7 +878,7 @@ def fit(self, X, feature_names):
 
         if len(base_term.penalties) > 1:
             raise NotImplementedError(
-                'bs="sz" currently requires a singly penalized base smooth.'
+                '"sz" smooth cannot use a multiply penalized basis (wrong basis in xt)'
             )
 
         self._base_term = base_term
diff --git a/nampy/gam/smooths/tensor/marginals.py b/nampy/gam/smooths/tensor/marginals.py
index dc11c7ce..04afb8fb 100644
--- a/nampy/gam/smooths/tensor/marginals.py
+++ b/nampy/gam/smooths/tensor/marginals.py
@@ -7,11 +7,12 @@
 from ...penalties.tensor import normalize_tensor_marginal_penalty
 from ..algebra import rowwise_kronecker
 from ..smooth_base import column_as_float
+from ..univariate.bs import DerivativeBSplineTerm1D
 from ..univariate.cr import CubicSplineTerm
 from ..univariate.ps import PSplineTerm1D
 from ..univariate.tp import ThinPlateSplineTerm
 
-TENSOR_MARGINAL_BASES = frozenset({"cr", "cs", "cc", "cp", "ps", "tp", "ts"})
+TENSOR_MARGINAL_BASES = frozenset({"bs", "cr", "cs", "cc", "cp", "ps", "tp", "ts"})
 
 
 def _as_marginal_features(feature):
@@ -93,6 +94,23 @@ def make_tensor_marginal_term(
             metadata=metadata,
         )
 
+    if basis == "bs":
+        if len(marginal_features) != 1:
+            raise ValueError("Tensor marginal basis 'bs' only handles one feature.")
+        return DerivativeBSplineTerm1D(
+            feature=marginal_features[0],
+            k=k,
+            m=m,
+            label=str(feature),
+            smoothing_id=None,
+            by=None,
+            select=False,
+            fixed=False,
+            constraint_mode=constraint_mode,
+            knots=knots,
+            metadata=metadata,
+        )
+
     if basis in {"tp", "ts"}:
         return ThinPlateSplineTerm(
             feature=marginal_features,
diff --git a/nampy/gam/smooths/univariate/__init__.py b/nampy/gam/smooths/univariate/__init__.py
index 3021ca13..98b81bf8 100644
--- a/nampy/gam/smooths/univariate/__init__.py
+++ b/nampy/gam/smooths/univariate/__init__.py
@@ -1,10 +1,17 @@
+from .bs import DerivativeBSplineTerm1D
 from .cr import CubicSplineTerm
 from .ps import PSplineTerm1D
 from .tp import ThinPlateSplineTerm
 
+bs = DerivativeBSplineTerm1D
 cr = cs = cc = CubicSplineTerm
 cp = ps = PSplineTerm1D
 tp = ts = ThinPlateSplineTerm
 
-__all__ = ["CubicSplineTerm", "PSplineTerm1D", "ThinPlateSplineTerm"]
-__all__ += ["cr", "cs", "cc", "cp", "ps", "tp", "ts"]
+__all__ = [
+    "DerivativeBSplineTerm1D",
+    "CubicSplineTerm",
+    "PSplineTerm1D",
+    "ThinPlateSplineTerm",
+]
+__all__ += ["bs", "cr", "cs", "cc", "cp", "ps", "tp", "ts"]
diff --git a/nampy/gam/smooths/univariate/bs.py b/nampy/gam/smooths/univariate/bs.py
new file mode 100644
index 00000000..b7a86c7d
--- /dev/null
+++ b/nampy/gam/smooths/univariate/bs.py
@@ -0,0 +1,352 @@
+"""Derivative-penalized B-spline smooth term (``bs='bs'``)."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from ...constraints.absorption import (
+    apply_linear_constraint,
+    should_apply_identifiability_constraint,
+)
+from ...penalties import (
+    PenaltySpec,
+    normalize_penalty_spec,
+    penalty_id_for_local_index,
+)
+from ...penalties.algebra import penalty_rescale_factor, scale_penalty
+from ...splines.univariate.bs import (
+    build_derivative_bspline_setup,
+    normalize_bspline_orders,
+    predict_derivative_bspline,
+)
+from ..registry import register_smooth
+from ..smooth_base import (
+    BaseSmoothTerm,
+    _resolve_feature,
+    build_penalty_definition,
+    by_values_from_new_data,
+    column_as_numeric_array,
+    linear_functional_basis,
+    linear_functional_by_state,
+)
+
+
+@register_smooth("bs")
+class DerivativeBSplineTerm1D(BaseSmoothTerm):
+    term_type = "smooth"
+    basis_name = "bs"
+    supports_tensor_marginal = False
+
+    def __init__(
+        self,
+        feature,
+        k=-1,
+        m=None,
+        label=None,
+        term_id=None,
+        smoothing_id=None,
+        by=None,
+        sp=None,
+        select=False,
+        fixed=False,
+        constraint_mode="auto",
+        pc=None,
+        knots=None,
+        null_penalty_tol=1e-10,
+        metadata=None,
+    ):
+        super().__init__(
+            feature=feature,
+            label=label,
+            term_id=term_id,
+            smoothing_id=smoothing_id,
+            by=by,
+            sp=sp,
+            metadata=metadata,
+        )
+        self.k = int(k)
+        self.m = normalize_bspline_orders(m)
+        self.select = bool(select)
+        self.fixed = bool(fixed)
+        self.constraint_mode = str(constraint_mode).lower()
+        self.pc = pc
+        self.knots = knots
+        self.null_penalty_tol = float(null_penalty_tol)
+
+        if self.select and self.fixed:
+            raise ValueError("select=True and fixed=True are incompatible.")
+        if self.constraint_mode not in {"auto", "factor_by", "always", "never"}:
+            raise ValueError(
+                "constraint_mode must be one of "
+                "{'auto', 'factor_by', 'always', 'never'}."
+            )
+
+        self._feature_index = None
+        self._feature_name = None
+        self._by_state = None
+        self._basis_train = None
+        self._penalties = None
+        self._setup = None
+        self._linear_functional = False
+
+    @property
+    def n_main_penalties(self):
+        return 0 if self.fixed else len(self.m) - 1
+
+    @property
+    def expected_linked_penalty_count(self):
+        return None if self.select else self.n_main_penalties
+
+    def _fit_constraint_policy(self, base, setup_base, penalties):
+        penalties_in = [] if self.fixed else list(penalties)
+        mode = self.constraint_mode
+        if mode == "factor_by" and not self._by_state.is_present:
+            raise ValueError(
+                "constraint_mode='factor_by' requires a numeric indicator `by` column."
+            )
+        should_constrain = (
+            mode == "factor_by"
+            or should_apply_identifiability_constraint(
+                self._by_state,
+                mode,
+                default_when_auto=True,
+            )
+        )
+        if should_constrain:
+            _, transformed, transform = apply_linear_constraint(
+                setup_base,
+                penalties_in,
+                np.asarray(setup_base, dtype=np.float64).mean(axis=0),
+            )
+            base_out = np.asarray(base, dtype=np.float64) @ transform
+            kind = "factor_by" if mode == "factor_by" else "sum_to_zero"
+        else:
+            transformed = penalties_in
+            transform = None
+            base_out = np.asarray(base, dtype=np.float64)
+            kind = None
+        base_out = self._apply_cached_by(base_out)
+        self._basis_train = np.asarray(base_out, dtype=np.float64)
+        self._penalties = [np.asarray(S, dtype=np.float64) for S in transformed]
+        self._record_constraint_result(
+            kind,
+            transform,
+            absorbed_by=("runtime" if transform is not None else None),
+        )
+
+    def fit(self, X, feature_names):
+        self._X_train = np.asarray(X, dtype=object).copy()
+        idx, feature_name = _resolve_feature(self.feature, feature_names)
+        self._feature_index = idx
+        self._feature_name = feature_name
+        self._set_resolved_features([feature_name])
+        x_values = column_as_numeric_array(X, idx)
+        self._set_by_state(X, feature_names)
+
+        self._linear_functional = np.asarray(x_values).ndim == 2
+        if self._linear_functional:
+            if self._by_state is None or np.asarray(self._by_state.values).ndim != 2:
+                raise ValueError(
+                    "B-spline linear-functional terms require matrix-valued by weights."
+                )
+            by_weights = np.asarray(self._by_state.values, dtype=np.float64)
+            if by_weights.shape != np.asarray(x_values).shape:
+                raise ValueError(
+                    "Linear-functional feature locations and by weights must have equal shape."
+                )
+            setup_values = np.asarray(x_values, dtype=np.float64).reshape(-1)
+            self._by_state = linear_functional_by_state(self._by_state)
+        else:
+            setup_values = np.asarray(x_values, dtype=np.float64).reshape(-1)
+
+        shared_X = self._linked_id_setup_matrix(feature_names)
+        pooled_setup = shared_X is not None
+        if pooled_setup:
+            if self._linear_functional:
+                raise NotImplementedError(
+                    "Linked-id pooling is not available for linear-functional B-splines."
+                )
+            setup_values = np.asarray(
+                column_as_numeric_array(shared_X, idx), dtype=np.float64
+            ).reshape(-1)
+
+        self._setup = build_derivative_bspline_setup(
+            setup_values,
+            feature_index=idx,
+            feature_name=feature_name,
+            bs_dim=self.k,
+            m=self.m,
+            knots=self.knots,
+        )
+        point_base = np.asarray(
+            predict_derivative_bspline(x_values, self._setup), dtype=np.float64
+        )
+        if self._linear_functional:
+            base = linear_functional_basis(
+                x_values,
+                by_weights,
+                lambda values: predict_derivative_bspline(values, self._setup),
+            )
+            setup_base = np.asarray(base, dtype=np.float64)
+        else:
+            base = point_base
+            setup_base = np.asarray(self._setup.basis_train, dtype=np.float64)
+
+        raw_penalties = [np.asarray(S, dtype=np.float64) for S in self._setup.penalties]
+        scales = [penalty_rescale_factor(setup_base, S) for S in raw_penalties]
+        self._set_penalty_rescale_factors(scales)
+        scaled_penalties = [scale_penalty(setup_base, S) for S in raw_penalties]
+
+        if self.pc is not None:
+            Bc, Sc, C, _ = self._apply_point_constraint(
+                base,
+                scaled_penalties,
+                self.pc,
+                feature_names=[self._feature_name],
+                point_basis_fn=lambda pts: predict_derivative_bspline(pts, self._setup)[
+                    0
+                ],
+                fixed=self.fixed,
+            )
+            self._basis_train = np.asarray(Bc, dtype=np.float64)
+            self._penalties = Sc
+            self._record_constraint_result("pc", C, absorbed_by="runtime")
+            return self
+
+        self._fit_constraint_policy(base, setup_base, scaled_penalties)
+        return self
+
+    @property
+    def basis_train(self):
+        self._require_fitted()
+        return self._basis_train
+
+    @property
+    def penalties(self):
+        self._require_fitted()
+        return self._penalties
+
+    @property
+    def n_coef(self):
+        self._require_fitted()
+        return int(self._basis_train.shape[1])
+
+    def get_penalty_definitions(self):
+        self._require_fitted()
+        raw = list(self.penalties)
+        if not raw:
+            return []
+        selection_defs = self._build_selection_penalty_definitions(
+            raw,
+            null_penalty_tol=self.null_penalty_tol,
+        )
+        sp_vals = self._normalized_term_sp(len(raw) + len(selection_defs))
+        definitions = []
+        for j, penalty in enumerate(raw):
+            sid = (
+                None
+                if self.smoothing_id is None
+                else penalty_id_for_local_index(
+                    self.smoothing_id, j, n_penalties=len(raw)
+                )
+            )
+            sp_j = sp_vals[j] if j < len(sp_vals) else None
+            definitions.append(
+                build_penalty_definition(
+                    self,
+                    penalty,
+                    kind="smooth",
+                    smoothing_id=sid,
+                    sp_value_in=sp_j,
+                    metadata_extra={
+                        "term_sp": sp_j,
+                        "m": self.m,
+                        "derivative_order": self.m[j + 1],
+                        "is_selection_penalty": False,
+                    },
+                    local_penalty_index=j,
+                )
+            )
+
+        for offset, selection in enumerate(selection_defs, start=len(raw)):
+            sp_j = sp_vals[offset] if offset < len(sp_vals) else None
+            if sp_j is None:
+                definitions.append(selection)
+                continue
+            definitions.append(
+                normalize_penalty_spec(
+                    PenaltySpec(
+                        matrix=np.asarray(selection.matrix, dtype=np.float64),
+                        smoothing_id=selection.smoothing_id,
+                        kind=selection.kind,
+                        rank=selection.rank,
+                        null_space_dim=selection.null_space_dim,
+                        is_null_space_penalty=selection.is_null_space_penalty,
+                        sp_mode="fixed" if sp_j >= 0 else "estimate",
+                        sp_value=float(sp_j) if sp_j >= 0 else None,
+                        metadata=dict(selection.metadata),
+                    )
+                )
+            )
+        return definitions
+
+    def transform_new(self, X_new):
+        self._require_fitted()
+        x_values = column_as_numeric_array(X_new, self._feature_index)
+        if self._linear_functional:
+            basis = linear_functional_basis(
+                x_values,
+                by_values_from_new_data(X_new, self._by_state),
+                lambda values: predict_derivative_bspline(values, self._setup),
+            )
+            if self.constraint_transform is not None:
+                basis = basis @ self.constraint_transform
+            return np.asarray(basis, dtype=np.float64)
+        basis = predict_derivative_bspline(x_values, self._setup)
+        return self._apply_constraint_transform_and_by(basis, X_new)
+
+    def derivative_matrix(self, X_new=None, order=1):
+        self._require_fitted()
+        order = int(order)
+        if order < 1 or order > int(self._setup.degree):
+            raise ValueError(f"order must be between 1 and {self._setup.degree}.")
+        if self._linear_functional:
+            raise NotImplementedError(
+                "Derivatives of linear-functional B-spline terms require an "
+                "explicit functional derivative and are not inferred."
+            )
+        source = self._X_train if X_new is None else X_new
+        x_values = column_as_numeric_array(source, self._feature_index)
+        basis = predict_derivative_bspline(x_values, self._setup, deriv=order)
+        return self._apply_constraint_transform_and_by(basis, source)
+
+    def tensor_marginal_fit_matrices(
+        self, *, centered=False, apply_np=False, x_train=None
+    ):
+        del apply_np, x_train
+        if len(self._setup.penalties) != 1:
+            raise NotImplementedError(
+                "Sorry, tensor products of smooths with multiple penalties are not "
+                "supported."
+            )
+        setup_base = np.asarray(self._setup.basis_train, dtype=np.float64)
+        raw_penalty = np.asarray(self._setup.penalties[0], dtype=np.float64)
+        if centered:
+            return super().tensor_marginal_fit_matrices(centered=True)
+        return setup_base, raw_penalty, None
+
+    def tensor_marginal_predict_matrix(
+        self, X_new, *, centered=False, np_transform=None
+    ):
+        if centered:
+            basis = np.asarray(self.transform_new(X_new), dtype=np.float64)
+        else:
+            x_values = column_as_numeric_array(X_new, self._feature_index)
+            if np.asarray(x_values).ndim != 1:
+                raise NotImplementedError(
+                    "Matrix-valued B-splines cannot be tensor marginals."
+                )
+            basis = predict_derivative_bspline(x_values, self._setup)
+        if np_transform is not None:
+            basis = basis @ np.asarray(np_transform, dtype=np.float64)
+        return np.asarray(basis, dtype=np.float64)
diff --git a/nampy/gam/specs/__init__.py b/nampy/gam/specs/__init__.py
index bf4fdd49..b4835990 100644
--- a/nampy/gam/specs/__init__.py
+++ b/nampy/gam/specs/__init__.py
@@ -8,6 +8,7 @@
     CubicRegressionSmoothSpec,
     CubicShrinkageSmoothSpec,
     CyclicCubicRegressionSmoothSpec,
+    DerivativeBSplineSmoothSpec,
     FactorSmoothInteractionSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
@@ -29,6 +30,7 @@
     "CubicRegressionSmoothSpec",
     "CubicShrinkageSmoothSpec",
     "CyclicCubicRegressionSmoothSpec",
+    "DerivativeBSplineSmoothSpec",
     "FactorSmoothInteractionSpec",
     "PSplineSmoothSpec",
     "RandomEffectSmoothSpec",
diff --git a/nampy/gam/specs/modeling.py b/nampy/gam/specs/modeling.py
index feee207f..0a77e08f 100644
--- a/nampy/gam/specs/modeling.py
+++ b/nampy/gam/specs/modeling.py
@@ -42,7 +42,7 @@ def make_predictor_specs(model, feature_names, *, knots=None):
                     metadata={},
                 )
             )
-        elif basis in {"ps", "cp"}:
+        elif basis in {"bs", "ps", "cp"}:
             main_terms.append(
                 TermSpec(
                     kind="smooth",
@@ -105,7 +105,7 @@ def make_predictor_specs(model, feature_names, *, knots=None):
         else:
             raise NotImplementedError(
                 "Automatic main-effect construction currently supports "
-                "{'cr','cs','cc','cp','ps','tp','ts','re'}, "
+                "{'bs','cr','cs','cc','cp','ps','tp','ts','re'}, "
                 f"got {model.basis!r}."
             )
 
diff --git a/nampy/gam/specs/smooth.py b/nampy/gam/specs/smooth.py
index 3c8d1b09..299ac872 100644
--- a/nampy/gam/specs/smooth.py
+++ b/nampy/gam/specs/smooth.py
@@ -50,6 +50,14 @@ class PSplineSmoothSpec(BaseSmoothSpec):
     pc: Any = None
 
 
+@dataclass(frozen=True)
+class DerivativeBSplineSmoothSpec(BaseSmoothSpec):
+    bs: str = "bs"
+    m: Any = None
+    constraint_mode: str = "auto"
+    pc: Any = None
+
+
 @dataclass(frozen=True)
 class ShapeConstrainedSmoothSpec(BaseSmoothSpec):
     """SCAM SCOP-spline specification for a named shape basis code."""
@@ -86,6 +94,7 @@ class RandomEffectSmoothSpec(BaseSmoothSpec):
 @dataclass(frozen=True)
 class FactorSmoothInteractionSpec(BaseSmoothSpec):
     bs: str = "fs"
+    m: Any = None
     xt: Any = None
     constraint_mode: str = "auto"
 
@@ -93,6 +102,7 @@ class FactorSmoothInteractionSpec(BaseSmoothSpec):
 @dataclass(frozen=True)
 class SumToZeroFactorSmoothSpec(BaseSmoothSpec):
     bs: str = "sz"
+    m: Any = None
     xt: Any = None
     constraint_mode: str = "auto"
 
@@ -121,6 +131,7 @@ class TensorInteractionSmoothSpec(BaseSmoothSpec):
 SmoothSpec = Union[
     CubicRegressionSmoothSpec,
     CyclicCubicRegressionSmoothSpec,
+    DerivativeBSplineSmoothSpec,
     CubicShrinkageSmoothSpec,
     PSplineSmoothSpec,
     ShapeConstrainedSmoothSpec,
diff --git a/nampy/gam/specs/smooth_build.py b/nampy/gam/specs/smooth_build.py
index d75339de..94a017a2 100644
--- a/nampy/gam/specs/smooth_build.py
+++ b/nampy/gam/specs/smooth_build.py
@@ -12,6 +12,7 @@
     CubicRegressionSmoothSpec,
     CubicShrinkageSmoothSpec,
     CyclicCubicRegressionSmoothSpec,
+    DerivativeBSplineSmoothSpec,
     FactorSmoothInteractionSpec,
     PSplineSmoothSpec,
     RandomEffectSmoothSpec,
@@ -98,6 +99,20 @@ def _build_s_ps(opts) -> PSplineSmoothSpec:
     )
 
 
+def _build_s_bs(opts) -> DerivativeBSplineSmoothSpec:
+    return DerivativeBSplineSmoothSpec(
+        special="s",
+        k=opts["k"],
+        fx=opts["fx"],
+        select=opts["select"],
+        sp=opts["sp"],
+        knots=opts["knots"],
+        m=opts["m"],
+        constraint_mode=opts["constraint_mode"],
+        pc=opts["pc"],
+    )
+
+
 def _build_s_shape(opts) -> ShapeConstrainedSmoothSpec:
     return ShapeConstrainedSmoothSpec(
         special="s",
@@ -162,6 +177,7 @@ def _build_s_fs(opts) -> FactorSmoothInteractionSpec:
         select=opts["select"],
         sp=opts["sp"],
         knots=opts["knots"],
+        m=opts["m"],
         xt=opts["xt"],
     )
 
@@ -174,6 +190,7 @@ def _build_s_sz(opts) -> SumToZeroFactorSmoothSpec:
         select=opts["select"],
         sp=opts["sp"],
         knots=opts["knots"],
+        m=opts["m"],
         xt=opts["xt"],
     )
 
@@ -184,6 +201,7 @@ def _build_s_sz(opts) -> SumToZeroFactorSmoothSpec:
     "cc": _build_s_cc,
     "ps": _build_s_ps,
     "cp": _build_s_ps,
+    "bs": _build_s_bs,
     "tp": _build_s_tp,
     "ts": _build_s_ts,
     "re": _build_s_re,
@@ -276,7 +294,7 @@ def _is_vector_fx(fx) -> bool:
     return fx is not None and not np.isscalar(fx)
 
 
-_PC_SUPPORTED_S_BASES = {"cc", "cp", "cr", "cs", "ps", "tp", "ts"}
+_PC_SUPPORTED_S_BASES = {"bs", "cc", "cp", "cr", "cs", "ps", "tp", "ts"}
 
 
 def _dispatch_smooth_spec_from_options(opts) -> SmoothSpec:
@@ -292,7 +310,7 @@ def _dispatch_smooth_spec_from_options(opts) -> SmoothSpec:
             raise NotImplementedError(
                 f"pc= is not supported for s(..., bs={merged['bs']!r}); "
                 "point constraints are only supported for bs in "
-                "{'cc', 'cp', 'cr', 'cs', 'ps', 'tp', 'ts'}."
+                "{'bs', 'cc', 'cp', 'cr', 'cs', 'ps', 'tp', 'ts'}."
             )
         return builder(merged)
     if has_pc and special_key not in {"te", "ti"}:
@@ -490,7 +508,7 @@ def _default_k_for_smooth(kind, basis, features, default_k):
         # mgcv::te()/ti() default k to 5^d per marginal. The current
         # Python tensor surface supports one feature per marginal, so d = 1.
         return [5] * len(features)
-    if str(basis).lower() in {"tp", "ts"}:
+    if str(basis).lower() in {"bs", "tp", "ts"}:
         # mgcv/R/smooth.r::s() leaves k = -1; smooth.construct.tp.smooth.spec
         # resolves the d-dependent default M + c(8, 27, 100)[min(d, 3)] at
         # construction time (mgcv/R/smooth.r:1316-1318). A flat default here
diff --git a/nampy/gam/splines/univariate/__init__.py b/nampy/gam/splines/univariate/__init__.py
index 65e53898..75026dc4 100644
--- a/nampy/gam/splines/univariate/__init__.py
+++ b/nampy/gam/splines/univariate/__init__.py
@@ -1,5 +1,14 @@
 """Canonical univariate spline setup surface."""
 
+from .bs import (
+    DerivativeBSplineSetup,
+    build_derivative_bspline_setup,
+    derivative_bspline_design,
+    derivative_bspline_knots,
+    derivative_penalty_root,
+    normalize_bspline_orders,
+    predict_derivative_bspline,
+)
 from .cr import (
     CubicSplines,
     add_full_rank_shrinkage,
@@ -24,6 +33,13 @@
 from .tp import build_tprs_term_setup, predict_tprs_term
 
 __all__ = [
+    "DerivativeBSplineSetup",
+    "build_derivative_bspline_setup",
+    "derivative_bspline_design",
+    "derivative_bspline_knots",
+    "derivative_penalty_root",
+    "normalize_bspline_orders",
+    "predict_derivative_bspline",
     "add_full_rank_shrinkage",
     "bspline_design_matrix",
     "cyclic_cubic_bd",
diff --git a/nampy/gam/splines/univariate/bs.py b/nampy/gam/splines/univariate/bs.py
new file mode 100644
index 00000000..c24950f4
--- /dev/null
+++ b/nampy/gam/splines/univariate/bs.py
@@ -0,0 +1,307 @@
+"""Derivative-penalized B-spline primitives for ``mgcv``'s ``bs='bs'``."""
+
+from __future__ import annotations
+
+import warnings
+from dataclasses import dataclass
+
+import numpy as np
+from scipy.linalg import cholesky_banded, solve
+
+from ...linalg import symmetrize_matrix
+from .ps import bspline_design_matrix, pspline_predict_matrix
+
+
+def _is_missing_order(value) -> bool:
+    if value is None:
+        return True
+    if isinstance(value, str) and value.strip().upper() == "NA":
+        return True
+    try:
+        return bool(np.asarray(value).ndim == 0 and np.isnan(float(value)))
+    except (TypeError, ValueError):
+        return False
+
+
+def _integer_order(value, *, position: int) -> int:
+    try:
+        numeric = float(value)
+    except (TypeError, ValueError) as exc:
+        raise ValueError(
+            "For bs='bs', m entries must be non-negative integers or NA."
+        ) from exc
+    if not np.isfinite(numeric) or numeric != np.rint(numeric) or numeric < 0:
+        raise ValueError("For bs='bs', m entries must be non-negative integers or NA.")
+    del position
+    return int(numeric)
+
+
+def normalize_bspline_orders(m) -> tuple[int, ...]:
+    """Normalize ``m`` exactly as ``smooth.construct.bs.smooth.spec``."""
+    if m is None:
+        return (3, 2)
+    if np.isscalar(m):
+        if _is_missing_order(m):
+            return (3, 2)
+        degree = _integer_order(m, position=0)
+        return (degree, max(0, degree - 1))
+
+    values = list(np.asarray(m, dtype=object).ravel())
+    if len(values) == 0:
+        raise ValueError("For bs='bs', m must contain a spline degree.")
+    if len(values) == 1:
+        return normalize_bspline_orders(values[0])
+
+    if _is_missing_order(values[0]):
+        if _is_missing_order(values[1]):
+            return (3, 2)
+        values[0] = _integer_order(values[1], position=1) + 1
+    if _is_missing_order(values[1]):
+        values[1] = max(0, _integer_order(values[0], position=0) - 1)
+
+    orders = tuple(_integer_order(value, position=i) for i, value in enumerate(values))
+    derivative_orders = orders[1:]
+    if len(set(derivative_orders)) < len(derivative_orders):
+        raise ValueError("multiple penalties of the same order is silly")
+    if any(order > orders[0] for order in derivative_orders):
+        raise ValueError("requested non-existent derivative in B-spline penalty")
+    return orders
+
+
+def derivative_bspline_knots(x, bs_dim, degree, supplied_knots=None):
+    """Port the automatic, endpoint, full, and four-knot constructor rules."""
+    values = np.asarray(x, dtype=np.float64).ravel()
+    degree = int(degree)
+    bs_dim = int(bs_dim)
+    nk = bs_dim - degree + 1
+    if nk <= 0:
+        raise ValueError("basis dimension too small for b-spline order")
+    expected = nk + 2 * degree
+
+    knots = None
+    if supplied_knots is not None:
+        knots = np.asarray(supplied_knots, dtype=np.float64).ravel()
+
+    if knots is not None and knots.size == 4 and knots.size < expected:
+        limits = np.sort(knots)
+        if nk <= 1:
+            raise ValueError(
+                "basis dimension too small for automatic knot construction"
+            )
+        dx = (limits[3] - limits[0]) / float(nk - 1)
+        lower_outer = limits[0] - dx * degree
+        upper_outer = limits[3] + dx * degree
+        lower = np.linspace(lower_outer, limits[0], degree + 1)
+        middle = (
+            np.linspace(limits[1], limits[2], max(0, nk - 2))
+            if nk > 2
+            else np.empty(0, dtype=np.float64)
+        )
+        upper = np.linspace(limits[3], upper_outer, degree + 1)
+        return np.concatenate([lower, middle, upper])
+
+    if knots is None or knots.size == 2:
+        if knots is None:
+            lower = float(np.min(values))
+            upper = float(np.max(values))
+        else:
+            lower = float(np.min(knots))
+            upper = float(np.max(knots))
+            if lower > np.min(values) or upper < np.max(values):
+                raise ValueError("knot range does not include data")
+        if nk <= 1:
+            raise ValueError(
+                "basis dimension too small for automatic knot construction"
+            )
+        width = upper - lower
+        lower -= width * 0.001
+        upper += width * 0.001
+        dx = (upper - lower) / float(nk - 1)
+        return np.linspace(
+            lower - dx * degree,
+            upper + dx * degree,
+            expected,
+        )
+
+    if knots.size != expected:
+        raise ValueError(f"there should be {expected} supplied knots")
+    if np.any(np.diff(knots) < 0):
+        raise ValueError("supplied bs knots must be nondecreasing")
+    return knots
+
+
+def derivative_bspline_design(x, knots, degree, deriv=0):
+    """Evaluate the constructor basis inside its effective knot interval."""
+    values = np.asarray(x, dtype=np.float64).ravel()
+    knots = np.asarray(knots, dtype=np.float64).ravel()
+    degree = int(degree)
+    deriv = int(deriv)
+    if deriv < 0 or deriv > degree:
+        raise ValueError("requested non-existent derivative in B-spline penalty")
+    lower = float(knots[degree])
+    upper = float(knots[knots.size - degree - 1])
+    if np.min(values) < lower or np.max(values) > upper:
+        raise ValueError("x out of range")
+    return bspline_design_matrix(
+        values,
+        knots,
+        degree=degree,
+        deriv=deriv,
+        extrapolate=True,
+    )
+
+
+def derivative_penalty_root(knots, degree, derivative_order):
+    """Port mgcv's exact band-Cholesky integrated-derivative penalty root."""
+    knots = np.asarray(knots, dtype=np.float64).ravel()
+    degree = int(degree)
+    derivative_order = int(derivative_order)
+    polynomial_degree = degree - derivative_order
+    if polynomial_degree < 0:
+        raise ValueError("requested non-existent derivative in B-spline penalty")
+
+    n_basis = int(knots.size - degree - 1)
+    interior = knots[degree : n_basis + 1]
+    widths = np.diff(interior)
+    if np.any(widths < 0):
+        raise ValueError("supplied bs knots must be nondecreasing")
+
+    if polynomial_degree == 0:
+        points = 0.5 * (interior[:-1] + interior[1:])
+        design = derivative_bspline_design(
+            points,
+            knots,
+            degree=degree,
+            deriv=derivative_order,
+        )
+        return np.sqrt(widths)[:, None] * design
+
+    steps = np.repeat(widths / polynomial_degree, polynomial_degree)
+    points = np.cumsum(np.concatenate(([interior[0]], steps)), dtype=np.float64)
+    points = np.clip(points, interior[0], interior[-1])
+    design = derivative_bspline_design(
+        points,
+        knots,
+        degree=degree,
+        deriv=derivative_order,
+    )
+
+    local_nodes = np.linspace(-1.0, 1.0, polynomial_degree + 1)
+    vandermonde = local_nodes[:, None] ** np.arange(polynomial_degree + 1)[None, :]
+    inverse_vandermonde = solve(
+        vandermonde,
+        np.eye(polynomial_degree + 1),
+        assume_a="gen",
+        check_finite=True,
+    )
+    powers = np.add.outer(
+        np.arange(polynomial_degree + 1),
+        np.arange(polynomial_degree + 1),
+    )
+    gram = np.where(powers % 2 == 0, 2.0 / (powers + 1.0), 0.0)
+    local_weight = inverse_vandermonde.T @ gram @ inverse_vandermonde
+
+    n_nodes = widths.size * polynomial_degree + 1
+    band = np.zeros((polynomial_degree + 1, n_nodes), dtype=np.float64)
+    for interval, width in enumerate(widths):
+        base = interval * polynomial_degree
+        scale = float(width) / 2.0
+        for column in range(polynomial_degree + 1):
+            for row in range(column, polynomial_degree + 1):
+                band[row - column, base + column] += scale * local_weight[column, row]
+
+    factor = cholesky_banded(
+        band,
+        lower=True,
+        overwrite_ab=False,
+        check_finite=True,
+    )
+    root = factor[0, :, None] * design
+    for offset in range(1, polynomial_degree + 1):
+        root[:-offset, :] += factor[offset, :-offset, None] * design[offset:, :]
+    return np.asarray(root, dtype=np.float64)
+
+
+@dataclass
+class DerivativeBSplineSetup:
+    feature_index: int
+    feature_name: str
+    degree: int
+    derivative_orders: tuple[int, ...]
+    knots: np.ndarray
+    basis_train: np.ndarray
+    penalty_roots: tuple[np.ndarray, ...]
+    penalties: tuple[np.ndarray, ...]
+    ranks: tuple[int, ...]
+    null_space_dim: int
+    bs_dim: int
+    orders: tuple[int, ...]
+
+
+def build_derivative_bspline_setup(
+    x,
+    *,
+    feature_index,
+    feature_name,
+    bs_dim,
+    m=None,
+    knots=None,
+):
+    """Build raw ``Bspline.smooth`` basis and integrated derivative penalties."""
+    values = np.asarray(x, dtype=np.float64).ravel()
+    orders = normalize_bspline_orders(m)
+    degree = int(orders[0])
+    derivative_orders = tuple(int(value) for value in orders[1:])
+    resolved_bs_dim = max(10, degree) if int(bs_dim) < 0 else int(bs_dim)
+    full_knots = derivative_bspline_knots(
+        values,
+        bs_dim=resolved_bs_dim,
+        degree=degree,
+        supplied_knots=knots,
+    )
+    basis = derivative_bspline_design(values, full_knots, degree, deriv=0)
+    if np.any(np.sum(basis, axis=0) == 0.0):
+        warnings.warn(
+            "there is *no* information about some basis coefficients",
+            stacklevel=2,
+        )
+    if np.unique(values).size < resolved_bs_dim:
+        warnings.warn(
+            "basis dimension is larger than number of unique covariates",
+            stacklevel=2,
+        )
+
+    roots = tuple(
+        derivative_penalty_root(full_knots, degree, order)
+        for order in derivative_orders
+    )
+    penalties = tuple(symmetrize_matrix(root.T @ root) for root in roots)
+    ranks = tuple(int(resolved_bs_dim - order) for order in derivative_orders)
+    null_space_dim = int(min(derivative_orders))
+    return DerivativeBSplineSetup(
+        feature_index=int(feature_index),
+        feature_name=str(feature_name),
+        degree=degree,
+        derivative_orders=derivative_orders,
+        knots=np.asarray(full_knots, dtype=np.float64),
+        basis_train=np.asarray(basis, dtype=np.float64),
+        penalty_roots=roots,
+        penalties=penalties,
+        ranks=ranks,
+        null_space_dim=null_space_dim,
+        bs_dim=int(resolved_bs_dim),
+        orders=orders,
+    )
+
+
+def predict_derivative_bspline(x_new, setup: DerivativeBSplineSetup, deriv=0):
+    """Match ``Predict.matrix.Bspline.smooth`` including linear extrapolation."""
+    return np.asarray(
+        pspline_predict_matrix(
+            x_new,
+            setup.knots,
+            basis_order=int(setup.degree) - 1,
+            deriv=int(deriv),
+        ),
+        dtype=np.float64,
+    )
diff --git a/pyproject.toml b/pyproject.toml
index 9b18482a..628bc0f9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -205,6 +205,7 @@ markers = [
     "smooth_cr: tests covering cubic regression smooths",
     "smooth_cs: tests covering shrinkage cubic regression smooths",
     "smooth_cc: tests covering cyclic cubic smooths",
+    "smooth_bs: tests covering integrated-derivative B-spline smooths",
     "smooth_cp: tests covering cyclic P-spline smooths",
     "smooth_ps: tests covering P-spline smooths",
     "smooth_tp: tests covering thin plate smooths",
diff --git a/tests/SUBSYSTEM_COVERAGE.md b/tests/SUBSYSTEM_COVERAGE.md
index 86b2501e..67225661 100644
--- a/tests/SUBSYSTEM_COVERAGE.md
+++ b/tests/SUBSYSTEM_COVERAGE.md
@@ -15,7 +15,7 @@ local development notes rather than duplicated here.
 | Subsystem | Primary owner(s) | Primary tests | Notes |
 | --- | --- | --- | --- |
 | Formula/spec parsing | `nampy/gam/formula/`, `nampy/gam/specs/` | `tests/parity/test_mgcv_formula_parse_parity.py` | Direct formula parity vs `mgcv`. |
-| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`) across wrapped prediction, tensor/factor-smooth combinations, and the upstream tensor-`m` wrong-length warning and zero fallback. |
+| Smooth constructors / raw basis owners | `nampy/gam/smooths/`, `nampy/gam/splines/` | `tests/smooths/test_mgcv_raw_constructor_parity.py`, `tests/smooths/test_mgcv_smoothcon_parity.py`, `tests/parity/test_gam_spec_build_owner_contracts.py`, `tests/parity/test_mgcv_cp_combinations_parity.py`, `tests/parity/test_mgcv_bs_combinations_parity.py` | Basis, penalty, constructor, and smoothCon surfaces, including cyclic P-splines (`cp`) and integrated-derivative B-splines (`bs`) across prediction, multi-penalty, tensor/factor-smooth combinations, and the upstream tensor-`m` wrong-length warning and zero fallback. |
 | `pc=` / linked `id=` routing | smooth metadata + linked basis owners | `tests/smooths/test_mgcv_pc_id_parity.py` | Localizes shared-smoothing and point-constraint issues, including `te`/`ti` point constraints. |
 | Design / pre-fit assembly | `nampy/gam/compiler/`, `nampy/gam/fit/penalized_system.py` | `tests/optimization/test_mgcv_gam_setup_assembly_parity.py`, `tests/optimization/test_mgcv_preoptimization_blocks_parity.py`, `tests/optimization/test_mgcv_preoptimization_reparam_parity.py` | Setup, blocks, and reparameterization parity, including one global shared-component block with overlapping linear-predictor indices. |
 | Term wrapping / by-variable / offset routing | predictor wrapping + compiled term owners | `tests/optimization/test_gam_term_wrapping_owner_contracts.py` | Localizes wrapped predictor blocks, offset routing, and general-family block ownership before broader prediction parity. |
diff --git a/tests/TAXONOMY.md b/tests/TAXONOMY.md
index 8694bedc..0c4e6aec 100644
--- a/tests/TAXONOMY.md
+++ b/tests/TAXONOMY.md
@@ -12,7 +12,7 @@ The GAM test suite is intentionally overlapping. The goal is fast subset runs an
 - `tests/`: shared helpers, marker inference, taxonomy registry, static reference fixtures, and parity-generation R scripts
 
 ## Taxonomy Axes
-- `smooth_<name>`: `cr`, `cs`, `cc`, `cp`, `ps`, `tp`, `ts`, `te`, `ti`, `fs`, `sz`, `re`
+- `smooth_<name>`: `bs`, `cr`, `cs`, `cc`, `cp`, `ps`, `tp`, `ts`, `te`, `ti`, `fs`, `sz`, `re`
 - `family_<name>`: `gaussian`, `binomial`, `poisson`, `gamma`, `negbin`, `gaulss`, `gammals`, `general`
 - `method_<name>`: `fixed`, `reml`, `ml`, `laml`, `gcv`, `ubre`
 - `link_<name>`: `identity`, `log`, `inverse`, `logit`, `probit`, `cloglog`, `cauchit`, `sqrt`
diff --git a/tests/_taxonomy_registry.py b/tests/_taxonomy_registry.py
index 0ca0c388..d949bb36 100644
--- a/tests/_taxonomy_registry.py
+++ b/tests/_taxonomy_registry.py
@@ -14,6 +14,7 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
 
 
 _SMOOTH_MARK_NAMES = {
+    "bs": "smooth_bs",
     "cr": "smooth_cr",
     "cs": "smooth_cs",
     "cc": "smooth_cc",
@@ -124,6 +125,11 @@ def _leaf(leaf_id: str, *nodeid_parts: str) -> LeafCoverageExpectation:
 }
 
 _PRIMARY_COVERAGE_BY_MARK = {
+    "smooth_bs": (
+        "tests/parity/test_mgcv_bs_combinations_parity.py",
+        "tests/smooths/test_mgcv_raw_constructor_parity.py",
+        "tests/smooths/test_mgcv_smoothcon_parity.py",
+    ),
     "smooth_cr": (
         "tests/parity/test_mgcv_snapshot_core_matrix.py",
         "tests/smooths/test_mgcv_smoothcon_parity.py",
diff --git a/tests/conftest.py b/tests/conftest.py
index ef31dc19..3bc32726 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -236,6 +236,10 @@ def _infer_marks_from_text(texts: list[str]) -> set[str]:
     for token, mark in _SMOOTH_MARK_NAMES.items():
         if token in {"te", "ti"}:
             continue
+        if token == "bs":
+            if f'bs="{token}"' in joined or f"bs='{token}'" in joined:
+                marks.add(mark)
+            continue
         if (
             f'bs="{token}"' in joined
             or f"bs='{token}'" in joined
diff --git a/tests/diagnostics/test_mgcv_k_check_parity.py b/tests/diagnostics/test_mgcv_k_check_parity.py
index 15d8d321..d1900f3d 100644
--- a/tests/diagnostics/test_mgcv_k_check_parity.py
+++ b/tests/diagnostics/test_mgcv_k_check_parity.py
@@ -269,6 +269,14 @@ class TestKCheckParity:
                 {"x0", "x1"},
                 1e-4,
             ),
+            (
+                lambda: _make_gaussian_data(seed=602, n=180),
+                'y ~ s(x0, bs="bs", k=8) + s(x1, bs="bs", k=8)',
+                "gaussian",
+                "REML",
+                {"x0", "x1"},
+                1e-4,
+            ),
             (
                 lambda: _make_gaussian_data(seed=123, n=180),
                 'y ~ te(x0, x1, bs=["cr","cr"], k=[5,5])',
@@ -283,6 +291,7 @@ class TestKCheckParity:
             "gaussian_cr_fixed",
             "gaussian_ps",
             "gaussian_cp",
+            "gaussian_bs",
             "gaussian_te",
         ],
     )
diff --git a/tests/mgcv_parity_utils.py b/tests/mgcv_parity_utils.py
index 48eb9ce8..fbd3c993 100644
--- a/tests/mgcv_parity_utils.py
+++ b/tests/mgcv_parity_utils.py
@@ -1608,6 +1608,11 @@ def _run_mgcv_raw_constructor(
       knots = pack_vector(sm$knots, "numeric"),
       m = pack_vector(sm$m, "integer")
     ),
+    "Bspline.smooth" = list(
+      knots = pack_vector(sm$knots, "numeric"),
+      m = pack_vector(sm$m, "numeric"),
+      D = if (is.null(sm$D)) list() else lapply(sm$D, pack_matrix)
+    ),
     "tprs.smooth" = list(
       Xu = pack_matrix(sm$Xu),
       UZ = pack_matrix(sm$UZ),
@@ -1857,12 +1862,11 @@ def _run_mgcv_smoothcon_predict_matrix(
     knots: dict | None = None,
     absorb_cons: bool = True,
     scale_penalty: bool = True,
+    deriv: int | None = None,
 ):
     smooth_expr_r = _normalize_python_formula_text(smooth_expr)
     knots_payload = _normalize_raw_constructor_knots(knots)
-    _cache_key = _mgcv_fixture_key(
-        "smoothcon_predict_matrix",
-        {
+    cache_parts = {
             "version": _SMOOTHCON_PREDICT_MATRIX_FIXTURE_VERSION,
             "data": _df_fixture_repr(data),
             "newdata": _df_fixture_repr(newdata),
@@ -1870,8 +1874,10 @@ def _run_mgcv_smoothcon_predict_matrix(
             "knots": json.dumps(knots_payload, sort_keys=True, default=str),
             "absorb_cons": bool(absorb_cons),
             "scale_penalty": bool(scale_penalty),
-        },
-    )
+    }
+    if deriv is not None:
+        cache_parts["deriv"] = int(deriv)
+    _cache_key = _mgcv_fixture_key("smoothcon_predict_matrix", cache_parts)
     cached = _mgcv_fixture_load(_cache_key)
     if cached is not None:
         return _decode_packed_matrix_payload(cached)
@@ -1892,8 +1898,8 @@ def _run_mgcv_smoothcon_predict_matrix(
 }
 out <- args[[3]]
 kn <- NULL
-if (length(args) >= 8 && nzchar(args[[8]])) {
-  kraw <- fromJSON(args[[8]], simplifyVector = FALSE)
+if (length(args) >= 7 && nzchar(args[[7]])) {
+  kraw <- fromJSON(args[[7]], simplifyVector = FALSE)
   kn <- lapply(kraw, function(v) {
     if (is.null(v)) return(NULL)
     vals <- unlist(v, recursive = TRUE, use.names = FALSE)
@@ -1923,6 +1929,7 @@ def _run_mgcv_smoothcon_predict_matrix(
   absorb.cons = absorb_cons,
   scale.penalty = scale_penalty
 )[[1]]
+if (length(args) >= 8 && nzchar(args[[8]])) sm$deriv <- as.integer(args[[8]])
 pm <- PredictMat(sm, newd)
 write_json(
   list(X = pack_matrix(pm)),
@@ -1957,6 +1964,7 @@ def _run_mgcv_smoothcon_predict_matrix(
                 "true" if absorb_cons else "false",
                 "true" if scale_penalty else "false",
                 knots_json,
+                "" if deriv is None else str(int(deriv)),
             ),
             check=True,
             cwd=_REPO_ROOT,
diff --git a/tests/parity/test_gam_pipeline_combination_matrix.py b/tests/parity/test_gam_pipeline_combination_matrix.py
index 3b2022b9..1df32e62 100644
--- a/tests/parity/test_gam_pipeline_combination_matrix.py
+++ b/tests/parity/test_gam_pipeline_combination_matrix.py
@@ -200,7 +200,7 @@ def test_stage_1_supported_numeric_factor_interactions_rebuild_on_newdata_like_m
     np.testing.assert_allclose(actual_se, np.asarray(expected["se"]).ravel(), atol=2e-8, rtol=2e-8)
 
 
-@pytest.mark.parametrize("basis", ["cr", "cs", "cc", "cp", "ps", "tp", "ts"])
+@pytest.mark.parametrize("basis", ["bs", "cr", "cs", "cc", "cp", "ps", "tp", "ts"])
 def test_stage_2_univariate_runtime_boundary_predictions_and_se_match_mgcv(basis):
     """Every supported univariate runtime matches behavior at and beyond fit bounds."""
     data = _formula_data(seed=903, n=90)[["y", "x0"]].rename(columns={"x0": "x"})
diff --git a/tests/parity/test_mgcv_bs_combinations_parity.py b/tests/parity/test_mgcv_bs_combinations_parity.py
new file mode 100644
index 00000000..97c4561d
--- /dev/null
+++ b/tests/parity/test_mgcv_bs_combinations_parity.py
@@ -0,0 +1,198 @@
+"""Integrated parity coverage for derivative-penalized B-splines (``bs='bs'``)."""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from nampy.gam import GAM
+from nampy.gam.splines.univariate.bs import build_derivative_bspline_setup
+from tests.mgcv_parity_utils import (
+    _assert_basic_mgcv_parity,
+    _fit_nampy_model,
+    _fit_nampy_snapshot,
+    _run_mgcv_snapshot,
+)
+
+
+def _bs_data(seed=251, n=190):
+    rng = np.random.default_rng(seed)
+    x0 = rng.uniform(-2.0, 2.0, size=n)
+    x1 = rng.uniform(-1.5, 2.5, size=n)
+    z = 0.8 + rng.uniform(-0.4, 0.7, size=n)
+    f = np.asarray(["a", "b", "c"], dtype=object)[np.arange(n) % 3]
+    f1 = np.asarray(["u", "v"], dtype=object)[np.arange(n) % 2]
+    y = (
+        0.2
+        + z * np.sin(1.2 * x0)
+        + 0.35 * x1**2
+        + 0.2 * (f == "b")
+        - 0.15 * (f1 == "v")
+        + rng.normal(scale=0.12, size=n)
+    )
+    return pd.DataFrame({"y": y, "x0": x0, "x1": x1, "z": z, "f": f, "f1": f1})
+
+
+def _assert_snapshot_fit(actual, expected, *, atol=4e-8):
+    for key in ("response", "link"):
+        np.testing.assert_allclose(
+            actual["predictions"][key],
+            expected["predictions"][key],
+            atol=atol,
+            rtol=atol,
+        )
+    np.testing.assert_allclose(
+        actual["fit"]["edf_total"], expected["fit"]["edf_total"], atol=atol, rtol=atol
+    )
+
+
+def test_bs_numeric_by_select_true_matches_mgcv():
+    data = _bs_data(seed=252)
+    formula = 'y ~ s(x0, by=z, bs="bs", k=9)'
+    actual = _fit_nampy_snapshot(data, formula, "gaussian", "REML", select=True)
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML", select=True)
+    assert len(actual["fit"]["smoothing_params"]) == 2
+    _assert_basic_mgcv_parity(
+        actual,
+        expected,
+        pred_atol=3e-6,
+        pred_rtol=3e-6,
+        sp_log_atol=5e-5,
+        criterion_atol=4e-8,
+    )
+
+
+def test_bs_factor_by_fixed_sp_matches_mgcv():
+    data = _bs_data(seed=253)
+    formula = 'y ~ s(x0, by=f, bs="bs", k=8, sp=0.7)'
+    actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed")
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+    _assert_snapshot_fit(actual, expected)
+
+
+def test_bs_linked_multi_penalty_terms_pool_basis_and_share_two_sp():
+    data = _bs_data(seed=254)
+    formula = (
+        'y ~ s(x0, bs="bs", k=9, m=[3,2,0], id="derivative", sp=[0.6,0.8])'
+        ' + s(x1, bs="bs", k=9, m=[3,2,0], id="derivative")'
+    )
+    model = _fit_nampy_model(data, formula, "gaussian", "fixed")
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+    assert len(model.smoothing_params) == 2
+    runtimes = [
+        term.predict_fn.__self__
+        for term in model.gam_result_.compiled_model.compiled_terms
+        if term.term_type == "smooth"
+    ]
+    assert len(runtimes) == 2
+    np.testing.assert_allclose(runtimes[0]._setup.knots, runtimes[1]._setup.knots)
+    pooled = np.concatenate([data["x0"].to_numpy(), data["x1"].to_numpy()])
+    assert runtimes[0]._setup.knots[3] <= np.min(pooled)
+    assert runtimes[0]._setup.knots[-4] >= np.max(pooled)
+    actual = model.parity_snapshot(X=data, include_covariances=True)
+    _assert_snapshot_fit(actual, expected, atol=8e-8)
+
+
+def test_bs_fixed_term_has_no_penalty_and_matches_mgcv():
+    data = _bs_data(seed=255)
+    formula = 'y ~ s(x0, bs="bs", k=8, fx=True)'
+    model = _fit_nampy_model(data, formula, "gaussian", "fixed")
+    assert model.gam_result_.compiled_model.compiled_penalties == ()
+    actual = model.parity_snapshot(X=data, include_covariances=True)
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+    _assert_snapshot_fit(actual, expected)
+
+
+@pytest.mark.parametrize(
+    "formula",
+    [
+        'y ~ te(x0, x1, bs=["bs","bs"], k=[5,6], m=[2,1], sp=[0.6,0.8])',
+        'y ~ ti(x0, x1, bs=["bs","bs"], k=[5,6], m=[2,1], sp=[0.6,0.8])',
+    ],
+    ids=["te", "ti"],
+)
+def test_bs_tensor_marginal_fixed_sp_fit_matches_mgcv(formula):
+    data = _bs_data(seed=256, n=170)
+    actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed")
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+    _assert_snapshot_fit(actual, expected, atol=3e-7)
+
+
+def test_bs_multiply_penalized_tensor_margin_is_rejected():
+    data = _bs_data(seed=257, n=80)
+    formula = 'y ~ te(x0, x1, bs=["bs","bs"], k=[6,6], m=[[3,2,1],[3,2]])'
+    with pytest.raises(NotImplementedError, match="multiple penalties"):
+        GAM(formula=formula).fit(data=data)
+
+
+@pytest.mark.parametrize(
+    "formula",
+    [
+        'y ~ s(f, x0, bs="fs", k=6, xt="bs", m=[3,2], sp=[0.7,0.9,1.1])',
+        'y ~ s(f, f1, x0, bs="sz", k=6, xt="bs", m=[3,2], id="shared", sp=0.7)',
+    ],
+    ids=["fs", "sz"],
+)
+def test_bs_factor_smooth_base_fixed_sp_fit_matches_mgcv(formula):
+    data = _bs_data(seed=258, n=180)
+    actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed")
+    expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+    _assert_snapshot_fit(actual, expected, atol=4e-7)
+
+
+@pytest.mark.parametrize("basis", ["fs", "sz"])
+def test_bs_multiply_penalized_factor_smooth_base_is_rejected(basis):
+    data = _bs_data(seed=259, n=80)
+    formula = f'y ~ s(f, x0, bs="{basis}", k=6, xt="bs", m=[3,2,1])'
+    with pytest.raises(NotImplementedError, match="multiply penalized basis"):
+        GAM(formula=formula).fit(data=data)
+
+
+def test_bs_array_api_and_persistence_preserve_extrapolated_prediction(tmp_path):
+    data = _bs_data(seed=260, n=150)
+    features = data[["x0", "x1"]]
+    model = GAM(
+        family="gaussian",
+        basis="bs",
+        k=8,
+        optimize_smoothing=False,
+        smoothing_params=[0.5, 0.8],
+    ).fit(X=features, y=data["y"].to_numpy(dtype=np.float64))
+    newdata = pd.DataFrame({"x0": [-4.0, 0.0, 4.0], "x1": [-3.0, 1.0, 5.0]})
+    expected = model.predict(newdata, type="link")
+    path = tmp_path / "bs.pkl"
+    model.save_model(path)
+    restored = GAM.load_model(path)
+    np.testing.assert_allclose(restored.predict(newdata, type="link"), expected)
+
+
+@pytest.mark.parametrize(
+    "formula,message",
+    [
+        ('y ~ s(x0, bs="bs", k=2, m=[3,2])', "basis dimension too small"),
+        ('y ~ s(x0, bs="bs", k=6, m=[3,4])', "non-existent derivative"),
+        ('y ~ s(x0, bs="bs", k=6, m=[3,2,2])', "multiple penalties"),
+        ('y ~ s(x0, bs="bs", k=6, m=[3.5,2])', "non-negative integers"),
+    ],
+)
+def test_bs_invalid_orders_fail_loudly(formula, message):
+    data = _bs_data(seed=261, n=40)
+    with pytest.raises(ValueError, match=message):
+        GAM(formula=formula).fit(data=data)
+
+
+def test_bs_knot_validation_and_unique_covariate_warning():
+    x = np.linspace(0.0, 1.0, 30)
+    kwargs = {
+        "feature_index": 0,
+        "feature_name": "x",
+        "bs_dim": 8,
+        "m": (3, 2),
+    }
+    with pytest.raises(ValueError, match="knot range does not include data"):
+        build_derivative_bspline_setup(x, knots=[0.2, 0.8], **kwargs)
+    with pytest.raises(ValueError, match="there should be 12 supplied knots"):
+        build_derivative_bspline_setup(x, knots=np.linspace(0.0, 1.0, 11), **kwargs)
+    with pytest.warns(UserWarning, match="larger than number of unique"):
+        build_derivative_bspline_setup(np.repeat([0.0, 0.5, 1.0], 10), **kwargs)
diff --git a/tests/reference_fixtures/mgcv/0eb0869399d4795146973e718109793bd5258f5ca7f71fc6131db8cba0f6a431.json.gz b/tests/reference_fixtures/mgcv/0eb0869399d4795146973e718109793bd5258f5ca7f71fc6131db8cba0f6a431.json.gz
new file mode 100644
index 00000000..5b3e0d0d
Binary files /dev/null and b/tests/reference_fixtures/mgcv/0eb0869399d4795146973e718109793bd5258f5ca7f71fc6131db8cba0f6a431.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/24d23117529abe4f1f653f198a36e2994cfda6e89ed2c43405575bd75ef77c1c.json.gz b/tests/reference_fixtures/mgcv/24d23117529abe4f1f653f198a36e2994cfda6e89ed2c43405575bd75ef77c1c.json.gz
new file mode 100644
index 00000000..2d4bfb5e
Binary files /dev/null and b/tests/reference_fixtures/mgcv/24d23117529abe4f1f653f198a36e2994cfda6e89ed2c43405575bd75ef77c1c.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/2a0a62f50d17c9661e8ca96d1fc9467367434e05566e16d54139ed80cff20e83.json.gz b/tests/reference_fixtures/mgcv/2a0a62f50d17c9661e8ca96d1fc9467367434e05566e16d54139ed80cff20e83.json.gz
new file mode 100644
index 00000000..8a480345
Binary files /dev/null and b/tests/reference_fixtures/mgcv/2a0a62f50d17c9661e8ca96d1fc9467367434e05566e16d54139ed80cff20e83.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/311fcd1aa6cc7c8cfc33d8f8242f977ac427276add537d3dd2b31c584badfe9f.json.gz b/tests/reference_fixtures/mgcv/311fcd1aa6cc7c8cfc33d8f8242f977ac427276add537d3dd2b31c584badfe9f.json.gz
new file mode 100644
index 00000000..37e1ced0
Binary files /dev/null and b/tests/reference_fixtures/mgcv/311fcd1aa6cc7c8cfc33d8f8242f977ac427276add537d3dd2b31c584badfe9f.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/34a0a3253aea4216ab0474ced37364186fb5d25689cab6af1882d95204615617.json.gz b/tests/reference_fixtures/mgcv/34a0a3253aea4216ab0474ced37364186fb5d25689cab6af1882d95204615617.json.gz
new file mode 100644
index 00000000..ca0c7ff1
Binary files /dev/null and b/tests/reference_fixtures/mgcv/34a0a3253aea4216ab0474ced37364186fb5d25689cab6af1882d95204615617.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/3a3ff1b50693da98494357265963c2a79df488f86d49164fbb135aef2733b21b.json.gz b/tests/reference_fixtures/mgcv/3a3ff1b50693da98494357265963c2a79df488f86d49164fbb135aef2733b21b.json.gz
new file mode 100644
index 00000000..7c6f7c50
Binary files /dev/null and b/tests/reference_fixtures/mgcv/3a3ff1b50693da98494357265963c2a79df488f86d49164fbb135aef2733b21b.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/404d26f43c31ab09fe57f7cea5afca9e64da6424aa2beda4ed20529d3732e8e2.json.gz b/tests/reference_fixtures/mgcv/404d26f43c31ab09fe57f7cea5afca9e64da6424aa2beda4ed20529d3732e8e2.json.gz
new file mode 100644
index 00000000..09396bf6
Binary files /dev/null and b/tests/reference_fixtures/mgcv/404d26f43c31ab09fe57f7cea5afca9e64da6424aa2beda4ed20529d3732e8e2.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/43dc98479029d05966e3929faa1533d5b942f212365ec76cd279fa25ff2f2227.json.gz b/tests/reference_fixtures/mgcv/43dc98479029d05966e3929faa1533d5b942f212365ec76cd279fa25ff2f2227.json.gz
new file mode 100644
index 00000000..bf483e2f
Binary files /dev/null and b/tests/reference_fixtures/mgcv/43dc98479029d05966e3929faa1533d5b942f212365ec76cd279fa25ff2f2227.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/46efa9439cec6ae56ac7f2a6cca4b6de9e8f3d05f2c5a3f858e6d6c33252a1da.json.gz b/tests/reference_fixtures/mgcv/46efa9439cec6ae56ac7f2a6cca4b6de9e8f3d05f2c5a3f858e6d6c33252a1da.json.gz
new file mode 100644
index 00000000..2e556bc6
Binary files /dev/null and b/tests/reference_fixtures/mgcv/46efa9439cec6ae56ac7f2a6cca4b6de9e8f3d05f2c5a3f858e6d6c33252a1da.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/49843d7820ec120e95551e9af717817d171cc61ca0bf4a3914da21f4b9d7372f.json.gz b/tests/reference_fixtures/mgcv/49843d7820ec120e95551e9af717817d171cc61ca0bf4a3914da21f4b9d7372f.json.gz
new file mode 100644
index 00000000..8c700975
Binary files /dev/null and b/tests/reference_fixtures/mgcv/49843d7820ec120e95551e9af717817d171cc61ca0bf4a3914da21f4b9d7372f.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/4d6da96dc0ee999e18f41776904261fea0590d81a415a1333b26cd195d5a9381.json.gz b/tests/reference_fixtures/mgcv/4d6da96dc0ee999e18f41776904261fea0590d81a415a1333b26cd195d5a9381.json.gz
new file mode 100644
index 00000000..9f0a9450
Binary files /dev/null and b/tests/reference_fixtures/mgcv/4d6da96dc0ee999e18f41776904261fea0590d81a415a1333b26cd195d5a9381.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/4fecae07a0c2dad33e46095fc23e9025d71092f8dbe49249d2eabe0b4f5107ae.json.gz b/tests/reference_fixtures/mgcv/4fecae07a0c2dad33e46095fc23e9025d71092f8dbe49249d2eabe0b4f5107ae.json.gz
new file mode 100644
index 00000000..90605a55
Binary files /dev/null and b/tests/reference_fixtures/mgcv/4fecae07a0c2dad33e46095fc23e9025d71092f8dbe49249d2eabe0b4f5107ae.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/594d420679a791fe6f4678f5737980f25c657527327d64d530cf791a50ac2cf0.json.gz b/tests/reference_fixtures/mgcv/594d420679a791fe6f4678f5737980f25c657527327d64d530cf791a50ac2cf0.json.gz
new file mode 100644
index 00000000..5800d0f3
Binary files /dev/null and b/tests/reference_fixtures/mgcv/594d420679a791fe6f4678f5737980f25c657527327d64d530cf791a50ac2cf0.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/6b23cfd13f5ce74918136f96c64c53b0738812b8688a106b63b1f21beb540207.json.gz b/tests/reference_fixtures/mgcv/6b23cfd13f5ce74918136f96c64c53b0738812b8688a106b63b1f21beb540207.json.gz
new file mode 100644
index 00000000..f66db789
Binary files /dev/null and b/tests/reference_fixtures/mgcv/6b23cfd13f5ce74918136f96c64c53b0738812b8688a106b63b1f21beb540207.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/6d33330e07ad46282bd60f02dc616e55a7fedbd81ba4219496f247fdd29e7869.json.gz b/tests/reference_fixtures/mgcv/6d33330e07ad46282bd60f02dc616e55a7fedbd81ba4219496f247fdd29e7869.json.gz
new file mode 100644
index 00000000..0ba01885
Binary files /dev/null and b/tests/reference_fixtures/mgcv/6d33330e07ad46282bd60f02dc616e55a7fedbd81ba4219496f247fdd29e7869.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/6eb7725af71b026d2355315b0dcd65f3447714ed52a389258dbb31570878da35.json.gz b/tests/reference_fixtures/mgcv/6eb7725af71b026d2355315b0dcd65f3447714ed52a389258dbb31570878da35.json.gz
new file mode 100644
index 00000000..f121c90f
Binary files /dev/null and b/tests/reference_fixtures/mgcv/6eb7725af71b026d2355315b0dcd65f3447714ed52a389258dbb31570878da35.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/6f9c852cdc1757bcdc2e1fe67c3986c38d8eb8a891fbd269261cc0bcd0eefae4.json.gz b/tests/reference_fixtures/mgcv/6f9c852cdc1757bcdc2e1fe67c3986c38d8eb8a891fbd269261cc0bcd0eefae4.json.gz
new file mode 100644
index 00000000..a4df7510
Binary files /dev/null and b/tests/reference_fixtures/mgcv/6f9c852cdc1757bcdc2e1fe67c3986c38d8eb8a891fbd269261cc0bcd0eefae4.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/7f9febf8c9729bd46b55a6f3676e69b0cfdfca08b97267d47aed3992f21c1126.json.gz b/tests/reference_fixtures/mgcv/7f9febf8c9729bd46b55a6f3676e69b0cfdfca08b97267d47aed3992f21c1126.json.gz
new file mode 100644
index 00000000..c6d63e7e
Binary files /dev/null and b/tests/reference_fixtures/mgcv/7f9febf8c9729bd46b55a6f3676e69b0cfdfca08b97267d47aed3992f21c1126.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/95172a7e8ab49ed6dbeca183145c5f2b5cd3a0f690876286f7628e42d465a5b1.json.gz b/tests/reference_fixtures/mgcv/95172a7e8ab49ed6dbeca183145c5f2b5cd3a0f690876286f7628e42d465a5b1.json.gz
new file mode 100644
index 00000000..dc9911ce
Binary files /dev/null and b/tests/reference_fixtures/mgcv/95172a7e8ab49ed6dbeca183145c5f2b5cd3a0f690876286f7628e42d465a5b1.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/9748170f4d19d3d80435b541d0931a6bc157f5e6ae60f1bf53bdbb375b1203e5.json.gz b/tests/reference_fixtures/mgcv/9748170f4d19d3d80435b541d0931a6bc157f5e6ae60f1bf53bdbb375b1203e5.json.gz
new file mode 100644
index 00000000..0db760ef
Binary files /dev/null and b/tests/reference_fixtures/mgcv/9748170f4d19d3d80435b541d0931a6bc157f5e6ae60f1bf53bdbb375b1203e5.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/98886e51bd9ed977725cb2cc886a728382085abc02b81ac5949fe38a4ce06452.json.gz b/tests/reference_fixtures/mgcv/98886e51bd9ed977725cb2cc886a728382085abc02b81ac5949fe38a4ce06452.json.gz
new file mode 100644
index 00000000..65bb6eaa
Binary files /dev/null and b/tests/reference_fixtures/mgcv/98886e51bd9ed977725cb2cc886a728382085abc02b81ac5949fe38a4ce06452.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/994d01254950bae174fac5ea504fe16c6b0be932600e381f9ca5bf4114684622.json.gz b/tests/reference_fixtures/mgcv/994d01254950bae174fac5ea504fe16c6b0be932600e381f9ca5bf4114684622.json.gz
new file mode 100644
index 00000000..b9fdee70
Binary files /dev/null and b/tests/reference_fixtures/mgcv/994d01254950bae174fac5ea504fe16c6b0be932600e381f9ca5bf4114684622.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/a35cb4a6d7c6e42d66464d57db90c86c46077dfcf32c4205b0de04af5f38d680.json.gz b/tests/reference_fixtures/mgcv/a35cb4a6d7c6e42d66464d57db90c86c46077dfcf32c4205b0de04af5f38d680.json.gz
new file mode 100644
index 00000000..1865b9c5
Binary files /dev/null and b/tests/reference_fixtures/mgcv/a35cb4a6d7c6e42d66464d57db90c86c46077dfcf32c4205b0de04af5f38d680.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/b1fd352d357da993e2c0363532a2b0c40cd83d2a278b35b7061241b8d8b84262.json.gz b/tests/reference_fixtures/mgcv/b1fd352d357da993e2c0363532a2b0c40cd83d2a278b35b7061241b8d8b84262.json.gz
new file mode 100644
index 00000000..56c427d0
Binary files /dev/null and b/tests/reference_fixtures/mgcv/b1fd352d357da993e2c0363532a2b0c40cd83d2a278b35b7061241b8d8b84262.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/b2ceb7c2a69f3a560d00b71cad32d70cdb0b4352d8780f54ebd2b5db463fd0e9.json.gz b/tests/reference_fixtures/mgcv/b2ceb7c2a69f3a560d00b71cad32d70cdb0b4352d8780f54ebd2b5db463fd0e9.json.gz
new file mode 100644
index 00000000..6f117352
Binary files /dev/null and b/tests/reference_fixtures/mgcv/b2ceb7c2a69f3a560d00b71cad32d70cdb0b4352d8780f54ebd2b5db463fd0e9.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/ba4d496ef9580e349889ec71c38cbab6e0310bf1c2471dd8cd41af9a3ea4a0bf.json.gz b/tests/reference_fixtures/mgcv/ba4d496ef9580e349889ec71c38cbab6e0310bf1c2471dd8cd41af9a3ea4a0bf.json.gz
new file mode 100644
index 00000000..69a9cdfc
Binary files /dev/null and b/tests/reference_fixtures/mgcv/ba4d496ef9580e349889ec71c38cbab6e0310bf1c2471dd8cd41af9a3ea4a0bf.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/bb025af6e41479ac8a4e7504e7be071dba464856d7b97e0aa7c5276c6e4a0979.json.gz b/tests/reference_fixtures/mgcv/bb025af6e41479ac8a4e7504e7be071dba464856d7b97e0aa7c5276c6e4a0979.json.gz
new file mode 100644
index 00000000..9a8982e2
Binary files /dev/null and b/tests/reference_fixtures/mgcv/bb025af6e41479ac8a4e7504e7be071dba464856d7b97e0aa7c5276c6e4a0979.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/e1e7fe072e7198c29b099f70d55e7387c1e9611bd93d5c5318cf38f9ff949384.json.gz b/tests/reference_fixtures/mgcv/e1e7fe072e7198c29b099f70d55e7387c1e9611bd93d5c5318cf38f9ff949384.json.gz
new file mode 100644
index 00000000..7580e810
Binary files /dev/null and b/tests/reference_fixtures/mgcv/e1e7fe072e7198c29b099f70d55e7387c1e9611bd93d5c5318cf38f9ff949384.json.gz differ
diff --git a/tests/reference_fixtures/mgcv/e49d64152f0c9a05bdc0ee1150db228809521dcf9424034e13a48afd921facb7.json.gz b/tests/reference_fixtures/mgcv/e49d64152f0c9a05bdc0ee1150db228809521dcf9424034e13a48afd921facb7.json.gz
new file mode 100644
index 00000000..22d8f832
Binary files /dev/null and b/tests/reference_fixtures/mgcv/e49d64152f0c9a05bdc0ee1150db228809521dcf9424034e13a48afd921facb7.json.gz differ
diff --git a/tests/scam/test_generic_transform_contracts.py b/tests/scam/test_generic_transform_contracts.py
index 4666ca0c..2cf127e6 100644
--- a/tests/scam/test_generic_transform_contracts.py
+++ b/tests/scam/test_generic_transform_contracts.py
@@ -137,7 +137,7 @@ def test_ar1_reml_rejects_missing_correlation_likelihood_terms():
         ).fit(data=data)
 
 
-@pytest.mark.parametrize("basis", ["ps", "cp"])
+@pytest.mark.parametrize("basis", ["bs", "ps", "cp"])
 def test_pspline_exposes_exact_derivative_provider_at_new_data(basis):
     x = np.linspace(-1.5, 2.0, 60)
     data = pd.DataFrame({"y": np.sin(x), "x": x})
@@ -165,7 +165,7 @@ def test_pspline_exposes_exact_derivative_provider_at_new_data(basis):
     assert derivative.derivative_matrix.shape[0] == len(new_data)
 
 
-@pytest.mark.parametrize("basis", ["ps", "cp", "cr", "cc"])
+@pytest.mark.parametrize("basis", ["bs", "ps", "cp", "cr", "cc"])
 def test_linear_functional_smooth_is_available_through_generic_by_contract(basis):
     rng = np.random.default_rng(91)
     locations = np.tile(np.linspace(-1.0, 1.0, 11), (28, 1))
diff --git a/tests/smooths/test_mgcv_raw_constructor_parity.py b/tests/smooths/test_mgcv_raw_constructor_parity.py
index 899ddc5b..b496ba83 100644
--- a/tests/smooths/test_mgcv_raw_constructor_parity.py
+++ b/tests/smooths/test_mgcv_raw_constructor_parity.py
@@ -26,6 +26,7 @@
 from nampy.gam.smooths.tensor.ti import (
     InteractionTensorProductSplineTerm,
 )
+from nampy.gam.smooths.univariate.bs import DerivativeBSplineTerm1D
 from nampy.gam.smooths.univariate.cr import CubicSplineTerm
 from nampy.gam.smooths.univariate.ps import PSplineTerm1D
 from nampy.gam.smooths.univariate.tp import ThinPlateSplineTerm
@@ -180,6 +181,15 @@ def _build(_data):
     return _build
 
 
+def _bspline_fixed_knots(column: str, values):
+    knots = np.asarray(values, dtype=np.float64)
+
+    def _build(_data):
+        return {str(column): knots.copy()}
+
+    return _build
+
+
 def _paired_feature_knots(columns, n_knots: int):
     cols = [str(col) for col in columns]
 
@@ -412,6 +422,64 @@ def _build_cp_case_matrix():
     ]
 
 
+def _build_bs_case_matrix():
+    return [
+        _case(
+            "bs_default_k_default_m",
+            _factory(_make_univariate_data, seed=141),
+            'y ~ s(x, bs="bs")',
+        ),
+        _case(
+            "bs_na_degree_inference",
+            _factory(_make_univariate_data, seed=142),
+            'y ~ s(x, bs="bs", k=8, m=c(NA, 1))',
+        ),
+        _case(
+            "bs_multi_derivative_penalties",
+            _factory(_make_univariate_data, seed=143),
+            'y ~ s(x, bs="bs", k=10, m=[3, 2, 1, 0])',
+        ),
+        _case(
+            "bs_two_limit_knots",
+            _factory(_make_univariate_data, seed=144),
+            'y ~ s(x, bs="bs", k=10, m=[3, 2])',
+            knots_factory=_bspline_fixed_knots("x", [-3.0, 3.0]),
+        ),
+        _case(
+            "bs_irregular_full_knots",
+            _factory(_make_univariate_data, seed=145),
+            'y ~ s(x, bs="bs", k=10, m=[3, 2, 0])',
+            knots_factory=_bspline_fixed_knots(
+                "x",
+                [
+                    -5.0,
+                    -4.0,
+                    -3.0,
+                    -2.5,
+                    -1.7,
+                    -1.0,
+                    -0.2,
+                    0.4,
+                    1.1,
+                    1.7,
+                    2.5,
+                    3.0,
+                    4.0,
+                    5.0,
+                ],
+            ),
+            atol=2e-9,
+        ),
+        _case(
+            "bs_special_four_knots",
+            _factory(_make_univariate_data, seed=146),
+            'y ~ s(x, bs="bs", k=10, m=[3, 1])',
+            knots_factory=_bspline_fixed_knots("x", [-3.0, -2.2, 2.2, 3.0]),
+            atol=2e-9,
+        ),
+    ]
+
+
 def _build_tprs_case_matrix():
     cases = []
     for basis, seed_base in [("tp", 60), ("ts", 80)]:
@@ -600,6 +668,7 @@ def _build_factor_smooth_case_matrix():
         ("cr", "cr"),
         ("cs", "cs"),
         ("cc", "cc"),
+        ("bs", "bs"),
         ("ps", {"bs": "ps", "m": 2, "k": 7}),
         ("cp", {"bs": "cp", "m": 2, "k": 7}),
         ("ts", "ts"),
@@ -610,7 +679,7 @@ def _build_factor_smooth_case_matrix():
                 f"fs_base_{label}",
                 _make_fs_data,
                 f'y ~ s(f, x, bs="fs", xt={repr(xt_spec)})',
-                atol=1e-8 if label in {"ps", "cp", "ts"} else 1e-10,
+                atol=1e-8 if label in {"bs", "ps", "cp", "ts"} else 1e-10,
             )
         )
         cases.append(
@@ -618,7 +687,7 @@ def _build_factor_smooth_case_matrix():
                 f"sz_base_{label}",
                 _make_sz_data,
                 f'y ~ s(f1, f2, x, bs="sz", k=6, xt={repr(xt_spec)})',
-                atol=1e-8 if label in {"ps", "cp", "ts"} else 1e-10,
+                atol=1e-8 if label in {"bs", "ps", "cp", "ts"} else 1e-10,
             )
         )
 
@@ -692,6 +761,18 @@ def _build_tensor_case_matrix():
                 'y ~ ti(x0, x1, bs=["cp", "cp"], k=[5, 6], m=[[2, 1], [2, 2]])',
                 atol=1e-8,
             ),
+            _case(
+                "te_bs_bs_m",
+                _factory(_make_gaussian_data, seed=820, n=90),
+                'y ~ te(x0, x1, bs=["bs", "bs"], k=[5, 6], m=[[3, 2], [2, 1]])',
+                atol=1e-8,
+            ),
+            _case(
+                "ti_bs_bs_m",
+                _factory(_make_gaussian_data, seed=821, n=90),
+                'y ~ ti(x0, x1, bs=["bs", "bs"], k=[5, 6], m=[[3, 2], [2, 1]])',
+                atol=1e-8,
+            ),
             _case(
                 "te_tp_ts_m",
                 _factory(_make_gaussian_data, seed=804, n=90),
@@ -756,6 +837,7 @@ def _build_tensor_case_matrix():
     *_build_cubic_case_matrix(),
     *_build_ps_case_matrix(),
     *_build_cp_case_matrix(),
+    *_build_bs_case_matrix(),
     *_build_tprs_case_matrix(),
     *_build_re_case_matrix(),
     *_build_factor_smooth_case_matrix(),
@@ -915,6 +997,22 @@ def _serialize_ps_raw(term):
     )
 
 
+def _serialize_bs_raw(term):
+    setup = term._setup
+    return _common_raw_state(
+        "Bspline.smooth",
+        np.asarray(setup.basis_train, dtype=np.float64),
+        [np.asarray(S, dtype=np.float64) for S in setup.penalties],
+        rank=_scalar_or_list(list(setup.ranks)),
+        null_space_dim=int(setup.null_space_dim),
+        extra={
+            "knots": np.asarray(setup.knots, dtype=np.float64),
+            "m": _scalar_or_list(list(setup.orders)),
+            "D": [np.asarray(D, dtype=np.float64) for D in setup.penalty_roots],
+        },
+    )
+
+
 def _serialize_tprs_raw(term):
     setup = term._setup
     B = np.asarray(setup.basis_train, dtype=np.float64)
@@ -1154,6 +1252,8 @@ def _serialize_ti_raw(term, X):
 
 
 def _serialize_term_raw(term, X):
+    if isinstance(term, DerivativeBSplineTerm1D):
+        return _serialize_bs_raw(term)
     if isinstance(term, CubicSplineTerm):
         return _serialize_cubic_raw(term, X)
     if isinstance(term, PSplineTerm1D):
diff --git a/tests/smooths/test_mgcv_smoothcon_parity.py b/tests/smooths/test_mgcv_smoothcon_parity.py
index 864a29be..54b15e01 100644
--- a/tests/smooths/test_mgcv_smoothcon_parity.py
+++ b/tests/smooths/test_mgcv_smoothcon_parity.py
@@ -15,6 +15,7 @@
 from nampy.gam.formula import extract_formula_terms, parse_gam_formula
 from nampy.gam.linalg import matrix_self_gram
 from nampy.gam.linalg import symmetric_spectrum as penalty_spectrum
+from nampy.gam.smooths.univariate.bs import DerivativeBSplineTerm1D
 from nampy.gam.smooths.univariate.cr import CubicSplineTerm
 from nampy.gam.specs.build import build_formula_model
 from tests._mgcv_snapshot_parity_shared import (
@@ -32,6 +33,7 @@
     _run_mgcv_smoothcon_matrix,
     _run_mgcv_smoothcon_matrix_unscaled,
     _run_mgcv_smoothcon_penalties,
+    _run_mgcv_smoothcon_predict_matrix,
     _run_mgcv_snapshot,
 )
 
@@ -862,6 +864,140 @@ def test_cp_reml_fit_matches_mgcv(self):
         )
 
 
+class TestDerivativeBSplineSmooth:
+    """Integrated-derivative B-spline (bs='bs') constructor and fit parity."""
+
+    @staticmethod
+    def _make_data(seed=191, n=190):
+        rng = np.random.default_rng(seed)
+        x = rng.uniform(-2.0, 2.0, size=n)
+        y = np.sin(1.3 * x) + 0.2 * x**2 + rng.normal(scale=0.12, size=n)
+        return pd.DataFrame({"y": y, "x": x})
+
+    def test_bs_smoothcon_basis_and_multiple_penalties_match_mgcv(self):
+        data = self._make_data()
+        formula = 'y ~ s(x, bs="bs", k=10, m=[3,2,1,0])'
+        smooth_expr = 's(x, bs="bs", k=10, m=c(3,2,1,0))'
+        design = _compile_formula_design(data, formula)
+        expected_basis = _run_mgcv_smoothcon_matrix(data, smooth_expr)
+        expected_penalties = _run_mgcv_smoothcon_penalties(
+            data, smooth_expr, absorb_cons=True, scale_penalty=True
+        )
+
+        np.testing.assert_allclose(
+            design.design_matrix, expected_basis["X"], atol=2e-10, rtol=2e-10
+        )
+        actual = [pb.matrix for pb in design.compiled_penalties]
+        assert len(actual) == len(expected_penalties["S"]) == 3
+        for actual_penalty, expected_penalty in zip(
+            actual, expected_penalties["S"], strict=True
+        ):
+            np.testing.assert_allclose(
+                actual_penalty, expected_penalty, atol=2e-10, rtol=2e-10
+            )
+
+    def test_bs_point_constraint_matches_mgcv(self):
+        data = self._make_data(seed=192)
+        formula = 'y ~ s(x, bs="bs", k=9, pc=0.0, sp=0.6)'
+        actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed")
+        expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+        np.testing.assert_allclose(
+            actual["predictions"]["response"],
+            expected["predictions"]["response"],
+            atol=2e-9,
+            rtol=2e-9,
+        )
+
+    def test_bs_multiple_fixed_sp_fit_matches_mgcv(self):
+        data = self._make_data(seed=193)
+        formula = 'y ~ s(x, bs="bs", k=10, m=[3,2,0], sp=[0.7,0.9])'
+        actual = _fit_nampy_snapshot(data, formula, "gaussian", "fixed")
+        expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+        np.testing.assert_allclose(
+            actual["predictions"]["response"],
+            expected["predictions"]["response"],
+            atol=2e-9,
+            rtol=2e-9,
+        )
+        np.testing.assert_allclose(
+            actual["fit"]["cov_bayes"],
+            expected["fit"]["cov_bayes"],
+            atol=2e-9,
+            rtol=2e-9,
+        )
+
+    def test_bs_reml_fit_matches_mgcv(self):
+        data = self._make_data(seed=194, n=210)
+        formula = 'y ~ s(x, bs="bs", k=11, m=[3,2])'
+        actual = _fit_nampy_snapshot(data, formula, "gaussian", "REML")
+        expected = _run_mgcv_snapshot(data, formula, "gaussian", "REML")
+        _assert_basic_mgcv_parity(
+            actual,
+            expected,
+            pred_atol=3e-8,
+            pred_rtol=3e-8,
+            sp_log_atol=3e-7,
+            criterion_atol=3e-8,
+        )
+
+    def test_bs_select_reml_matches_mgcv(self):
+        data = self._make_data(seed=195, n=210)
+        formula = 'y ~ s(x, bs="bs", k=10)'
+        actual = _fit_nampy_snapshot(
+            data, formula, "gaussian", "REML", select=True
+        )
+        expected = _run_mgcv_snapshot(
+            data, formula, "gaussian", "REML", select=True
+        )
+        assert len(actual["fit"]["smoothing_params"]) == 2
+        _assert_basic_mgcv_parity(
+            actual,
+            expected,
+            pred_atol=3e-8,
+            pred_rtol=3e-8,
+            sp_log_atol=4e-7,
+            criterion_atol=3e-8,
+        )
+
+    def test_bs_derivative_prediction_and_linear_tails_match_mgcv(self):
+        data = self._make_data(seed=196)
+        newdata = pd.DataFrame({"x": [-3.0, -1.0, 0.5, 3.0]})
+        term = DerivativeBSplineTerm1D(feature="x", k=10, m=(3, 2))
+        term.fit(data[["x"]].to_numpy(dtype=np.float64), ["x"])
+        for order in (1, 2):
+            actual = term.derivative_matrix(
+                newdata[["x"]].to_numpy(dtype=np.float64), order=order
+            )
+            expected = _run_mgcv_smoothcon_predict_matrix(
+                data,
+                newdata,
+                's(x, bs="bs", k=10, m=c(3,2))',
+                deriv=order,
+            )
+            np.testing.assert_allclose(
+                actual, expected["X"], atol=2e-10, rtol=2e-10
+            )
+
+    def test_bs_four_knot_prediction_interval_matches_mgcv(self):
+        data = self._make_data(seed=197)
+        knots = {"x": [-3.0, -2.2, 2.2, 3.0]}
+        newdata = pd.DataFrame({"x": [-4.0, -2.7, 0.0, 2.7, 4.0]})
+        term = DerivativeBSplineTerm1D(
+            feature="x", k=10, m=(3, 1), knots=knots["x"]
+        )
+        term.fit(data[["x"]].to_numpy(dtype=np.float64), ["x"])
+        actual = term.transform_new(newdata[["x"]].to_numpy(dtype=np.float64))
+        expected = _run_mgcv_smoothcon_predict_matrix(
+            data,
+            newdata,
+            's(x, bs="bs", k=10, m=c(3,1))',
+            knots=knots,
+        )
+        np.testing.assert_allclose(
+            actual, expected["X"], atol=2e-10, rtol=2e-10
+        )
+
+
 class TestPSplineSmooth(_SharedTestPSplineSmooth):
     """P-spline (bs='ps') standalone parity against mgcv."""