diff --git a/README.md b/README.md
index c5412d0..7f05e2f 100644
--- a/README.md
+++ b/README.md
@@ -17,6 +17,7 @@ stay pure-Python while still offering compiled, numba-free fast paths.
 | Estimator | Compute API | Notes |
 |-----------|-------------|-------|
 | MiniRocket (univariate and multivariate) | `rocket.rocket_fit` / `rocket_transform` | numba-free; 2D or 3D NumPy input; equivalent to `MiniRocket` / `MiniRocketMultivariate`, no JIT warmup |
+| Multivariate MultiRocket | `sktime_cython.transformations.rocket.multirocket_fit` / `multirocket_transform` | numba-free; equivalent to `MultiRocketMultivariate`, no JIT warmup |
 
 Logic for further `sktime` cython based estimators should be added in this package,
 for a recipe, see the [extension guide](https://github.com/sktime/sktime-cython#adding-an-estimator) below.
diff --git a/setup.py b/setup.py
index ca61439..cc3fcab 100644
--- a/setup.py
+++ b/setup.py
@@ -12,16 +12,22 @@
 else:
     _fast = ["-O3", "-ffast-math"]
 
+# Fully qualified names, so extensions can live in any sktime_cython subpackage;
+# each source is the .pyx at the module's dotted path.
+_MODULES = [
+    "sktime_cython.transformations.rocket._minirocket_multivariate_cython",
+    "sktime_cython.transformations.rocket._multirocket_multivariate_cython",
+]
+
 extensions = [
     Extension(
-        "sktime_cython.transformations.rocket._minirocket_multivariate_cython",
-        sources=[
-            "sktime_cython/transformations/rocket/_minirocket_multivariate_cython.pyx"
-        ],
+        name,
+        sources=[name.replace(".", "/") + ".pyx"],
         include_dirs=[np.get_include()],
         define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")],
         extra_compile_args=_fast,
-    ),
+    )
+    for name in _MODULES
 ]
 
 setup(
diff --git a/sktime_cython/transformations/rocket/__init__.py b/sktime_cython/transformations/rocket/__init__.py
index 4a54dfc..fc86b03 100644
--- a/sktime_cython/transformations/rocket/__init__.py
+++ b/sktime_cython/transformations/rocket/__init__.py
@@ -4,5 +4,14 @@
     rocket_fit,
     rocket_transform,
 )
+from sktime_cython.transformations.rocket._multirocket import (
+    multirocket_fit,
+    multirocket_transform,
+)
 
-__all__ = ["rocket_fit", "rocket_transform"]
+__all__ = [
+    "multirocket_fit",
+    "multirocket_transform",
+    "rocket_fit",
+    "rocket_transform",
+]
diff --git a/sktime_cython/transformations/rocket/_common.py b/sktime_cython/transformations/rocket/_common.py
new file mode 100644
index 0000000..7a62517
--- /dev/null
+++ b/sktime_cython/transformations/rocket/_common.py
@@ -0,0 +1,116 @@
+"""Scaffolding shared by the MiniRocket and MultiRocket transforms.
+
+The dilation schedule and quantile points are identical for both transforms;
+the per-transform feature accounting happens in the callers. Pure numpy: these
+run once per fit on small arrays, so there is nothing here worth compiling.
+"""
+
+import numpy as np
+
+_NUM_KERNELS = 84
+
+
+def _fit_dilations(n_timepoints, num_features, max_dilations_per_kernel):
+    """Dilation schedule (pure numpy, copied from sktime numba module)."""
+    num_kernels = _NUM_KERNELS
+    if num_features < num_kernels:
+        num_features = num_kernels
+
+    num_features_per_kernel = num_features // num_kernels
+    true_max_dilations_per_kernel = min(
+        num_features_per_kernel, max_dilations_per_kernel
+    )
+    multiplier = num_features_per_kernel / true_max_dilations_per_kernel
+
+    max_exponent = np.log2((n_timepoints - 1) / (9 - 1))
+    dilations, num_features_per_dilation = np.unique(
+        np.logspace(0, max_exponent, true_max_dilations_per_kernel, base=2).astype(
+            np.int32
+        ),
+        return_counts=True,
+    )
+    num_features_per_dilation = (num_features_per_dilation * multiplier).astype(
+        np.int32
+    )
+
+    remainder = num_features_per_kernel - np.sum(num_features_per_dilation)
+    i = 0
+    while remainder > 0:
+        num_features_per_dilation[i] += 1
+        remainder -= 1
+        i = (i + 1) % len(num_features_per_dilation)
+
+    return dilations, num_features_per_dilation
+
+
+def _quantiles(n):
+    """Evenly-spaced low-discrepancy quantile points (copied from sktime)."""
+    return np.array(
+        [(_ * ((np.sqrt(5) + 1) / 2)) % 1 for _ in range(1, n + 1)], dtype=np.float32
+    )
+
+
+def _biases_from_C(C, quantiles, num_features_per_dilation, num_kernels):
+    """Take the per-combination quantiles of C, vectorized.
+
+    Equivalent to looping over combinations in dilation-major order and
+    calling np.quantile on each row with its own slice of quantiles.
+
+    Parameters
+    ----------
+    C : per-combination convolution outputs, one row per combination
+    quantiles : one quantile per output feature, in feature order
+    num_features_per_dilation : number of features per dilation
+    num_kernels : number of kernels per dilation
+
+    Returns
+    -------
+    the bias per output feature
+    """
+    features_per_combination = np.repeat(num_features_per_dilation, num_kernels)
+    combination = np.repeat(
+        np.arange(len(features_per_combination)), features_per_combination
+    )
+    # numba's np.quantile, op for op (numba/np/arraymath.py): it promotes both
+    # operands to float64, scales the quantile to a percentile and back, and
+    # interpolates between ranks f-1 and f. Matching it exactly matters because
+    # a bias landing on a repeated convolution value would otherwise
+    # reclassify a whole tied group of samples.
+    S = np.sort(C, axis=1).astype(np.float64)
+    n_timepoints = S.shape[1]
+    percentile = quantiles.astype(np.float64) * 100.0
+    rank = 1.0 + (n_timepoints - 1) * (percentile / 100.0)
+    floor = np.floor(rank)
+    weight = rank - floor
+    lower = (floor - 1.0).astype(np.intp)
+    upper = np.minimum(lower + 1, n_timepoints - 1)
+    return (
+        S[combination, lower] * (1 - weight) + S[combination, upper] * weight
+    ).astype(np.float32)
+
+
+def _check_dilations(dilations, n_timepoints, name):
+    """Check the in-bounds invariants the kernel relies on.
+
+    The Cython kernel runs with bounds checking off and integer division
+    unguarded, so a hand-built parameter tuple could otherwise corrupt memory
+    or divide by zero. ``_fit_dilations`` guarantees both conditions.
+
+    Parameters
+    ----------
+    dilations : the dilations for one representation
+    n_timepoints : length of the series that representation is applied to
+    name : parameter name, for the error message
+
+    Raises
+    ------
+    ValueError
+        if a dilation is below 1, or too large for the series length
+    """
+    if dilations.min() < 1:
+        raise ValueError(f"{name} must be >= 1, but found {dilations.min()}")
+    if 8 * dilations.max() >= n_timepoints:
+        raise ValueError(
+            f"{name} must satisfy 8 * dilation < n_timepoints={n_timepoints}, "
+            f"but found {dilations.max()}"
+        )
diff --git a/sktime_cython/transformations/rocket/_minirocket.py b/sktime_cython/transformations/rocket/_minirocket.py
index d0bf222..cafaa50 100644
--- a/sktime_cython/transformations/rocket/_minirocket.py
+++ b/sktime_cython/transformations/rocket/_minirocket.py
@@ -15,51 +15,16 @@
 import numpy as np
 
 from sktime_cython.transformations.rocket import _minirocket_multivariate_cython as _cy
+from sktime_cython.transformations.rocket._common import (
+    _NUM_KERNELS,
+    _biases_from_C,
+    _check_dilations,
+    _fit_dilations,
+    _quantiles,
+)
 
 __all__ = ["rocket_fit", "rocket_transform"]
 
-_NUM_KERNELS = 84
-
-
-def _fit_dilations(n_timepoints, num_features, max_dilations_per_kernel):
-    """Dilation schedule (pure numpy, copied from sktime numba module)."""
-    num_kernels = _NUM_KERNELS
-    if num_features < num_kernels:
-        num_features = num_kernels
-
-    num_features_per_kernel = num_features // num_kernels
-    true_max_dilations_per_kernel = min(
-        num_features_per_kernel, max_dilations_per_kernel
-    )
-    multiplier = num_features_per_kernel / true_max_dilations_per_kernel
-
-    max_exponent = np.log2((n_timepoints - 1) / (9 - 1))
-    dilations, num_features_per_dilation = np.unique(
-        np.logspace(0, max_exponent, true_max_dilations_per_kernel, base=2).astype(
-            np.int32
-        ),
-        return_counts=True,
-    )
-    num_features_per_dilation = (num_features_per_dilation * multiplier).astype(
-        np.int32
-    )
-
-    remainder = num_features_per_kernel - np.sum(num_features_per_dilation)
-    i = 0
-    while remainder > 0:
-        num_features_per_dilation[i] += 1
-        remainder -= 1
-        i = (i + 1) % len(num_features_per_dilation)
-
-    return dilations, num_features_per_dilation
-
-
-def _quantiles(n):
-    """Evenly-spaced low-discrepancy quantile points (copied from sktime)."""
-    return np.array(
-        [(_ * ((np.sqrt(5) + 1) / 2)) % 1 for _ in range(1, n + 1)], dtype=np.float32
-    )
-
 
 def _as_panel(X):
     """Return X as contiguous 3D float32/float64; 2D is a univariate panel."""
@@ -157,19 +122,7 @@ def rocket_fit(X, num_kernels=10_000, max_dilations_per_kernel=32, random_state=
         instance_indices,
     )
 
-    biases = np.zeros(num_kernels_ * int(num_features_per_kernel), dtype=np.float32)
-    feature_index_start = 0
-    combination_index = 0
-    for dilation_index in range(num_dilations):
-        nfd = num_features_per_dilation[dilation_index]
-        for _kernel_index in range(num_kernels_):
-            feature_index_end = feature_index_start + nfd
-            biases[feature_index_start:feature_index_end] = np.quantile(
-                C[combination_index],
-                quantiles[feature_index_start:feature_index_end],
-            ).astype(np.float32)
-            feature_index_start = feature_index_end
-            combination_index += 1
+    biases = _biases_from_C(C, quantiles, num_features_per_dilation, num_kernels_)
 
     return (
         num_channels_per_combination,
@@ -200,6 +153,7 @@ def rocket_transform(X, parameters, n_jobs=1):
     np.ndarray, shape (n_instances, n_features), float32
     """
     X = _as_panel(X)
+    _check_dilations(parameters[2], X.shape[2], "dilations")
     n_instances = X.shape[0]
 
     if n_jobs < 1 or n_jobs > multiprocessing.cpu_count():
diff --git a/sktime_cython/transformations/rocket/_multirocket.py b/sktime_cython/transformations/rocket/_multirocket.py
new file mode 100644
index 0000000..92a72bf
--- /dev/null
+++ b/sktime_cython/transformations/rocket/_multirocket.py
@@ -0,0 +1,221 @@
+"""Numba-free MultiRocket multivariate transform (pure-numpy + Cython core).
+
+Compute layer with no sktime dependency: numpy arrays in, numpy arrays out.
+``fit`` returns a parameter tuple; ``transform`` applies it. An sktime
+``BaseTransformer`` wrapper delegates to these two functions.
+
+MultiRocket applies the MiniRocket kernels to both the raw series and its
+first-order difference, and takes four features per kernel rather than one.
+"""
+
+import multiprocessing
+from concurrent.futures import ThreadPoolExecutor
+
+import numpy as np
+
+from sktime_cython.transformations.rocket import _multirocket_multivariate_cython as _cy
+from sktime_cython.transformations.rocket._common import (
+    _NUM_KERNELS,
+    _biases_from_C,
+    _check_dilations,
+    _fit_dilations,
+    _quantiles,
+)
+
+__all__ = ["multirocket_fit", "multirocket_transform"]
+
+
+def _normalise(X):
+    """Normalise each series to zero mean and unit variance."""
+    return (X - X.mean(axis=-1, keepdims=True)) / (X.std(axis=-1, keepdims=True) + 1e-8)
+
+
+def _get_parameter(X, num_kernels, max_dilations_per_kernel, seed):
+    """Fit dilations, channel selections, and biases for one representation."""
+    if seed is not None:
+        np.random.seed(seed)
+
+    n_instances, n_columns, n_timepoints = X.shape
+    num_kernels_ = _NUM_KERNELS
+
+    dilations, num_features_per_dilation = _fit_dilations(
+        n_timepoints, num_kernels, max_dilations_per_kernel
+    )
+    num_features_per_kernel = np.sum(num_features_per_dilation)
+    quantiles = _quantiles(num_kernels_ * num_features_per_kernel)
+
+    num_combinations = num_kernels_ * len(dilations)
+
+    max_num_channels = min(n_columns, 9)
+    max_exponent = np.log2(max_num_channels + 1)
+
+    num_channels_per_combination = (
+        2 ** np.random.uniform(0, max_exponent, num_combinations)
+    ).astype(np.int32)
+
+    channel_indices = np.zeros(num_channels_per_combination.sum(), dtype=np.int32)
+    num_channels_start = 0
+    for combination_index in range(num_combinations):
+        n_this = num_channels_per_combination[combination_index]
+        num_channels_end = num_channels_start + n_this
+        channel_indices[num_channels_start:num_channels_end] = np.random.choice(
+            n_columns, n_this, replace=False
+        )
+        num_channels_start = num_channels_end
+
+    # biases: re-seed (matching numba _fit_biases), draw one instance index per
+    # combination, build C in Cython, then quantile per combination.
+    if seed is not None:
+        np.random.seed(seed)
+    instance_indices = np.array(
+        [np.random.randint(n_instances) for _ in range(num_combinations)],
+        dtype=np.int32,
+    )
+    C = _cy.fit_biases(
+        X,
+        num_channels_per_combination,
+        channel_indices,
+        dilations.astype(np.int32),
+        num_features_per_dilation.astype(np.int32),
+        instance_indices,
+    )
+    biases = _biases_from_C(C, quantiles, num_features_per_dilation, num_kernels_)
+
+    return (
+        num_channels_per_combination,
+        channel_indices,
+        dilations.astype(np.int32),
+        num_features_per_dilation.astype(np.int32),
+        biases,
+    )
+
+
+def multirocket_fit(
+    X,
+    num_kernels=6_250,
+    max_dilations_per_kernel=32,
+    normalise=False,
+    random_state=None,
+):
+    """Fit dilations, channel selections, and biases, on X and its difference.
+
+    Parameters
+    ----------
+    X : 3D np.ndarray, shape (n_instances, n_columns, n_timepoints)
+        panel of time series; normalised in its own dtype, then cast to float64
+        (as in sktime).
+    num_kernels : int, default=6250
+        number of kernels; rounded down to a multiple of 84 (min 84).
+    max_dilations_per_kernel : int, default=32
+        maximum number of dilations per kernel.
+    normalise : bool, default=False
+        whether to normalise each series before fitting.
+    random_state : int or None, default=None
+        seed for reproducibility.
+
+    Returns
+    -------
+    parameters : flattened tuple of tuple of np.ndarray
+        the parameters for the raw and the differenced representation, ready
+        to pass to ``transform``.
+    """
+    if random_state is not None and not isinstance(random_state, (int, np.integer)):
+        raise ValueError(
+            f"random_state must be int or None, but found {type(random_state)}"
+        )
+    seed = np.int32(random_state) if random_state is not None else None
+
+    # sktime's order: normalise in the input dtype, pad, then cast to float64
+    if normalise:
+        X = _normalise(X)
+
+    if X.shape[2] < 10:
+        # handling very short series (like PensDigit from the MTSC archive)
+        # series have to be at least a length of 10 (including differencing)
+        padded = np.zeros((X.shape[0], X.shape[1], 10), dtype=X.dtype)
+        padded[:, :, : X.shape[2]] = X
+        X = padded
+    X = np.ascontiguousarray(X, dtype=np.float64)
+
+    args = (num_kernels, max_dilations_per_kernel, seed)
+    return (
+        _get_parameter(X, *args),
+        _get_parameter(np.ascontiguousarray(np.diff(X, 1)), *args),
+    )
+
+
+def multirocket_transform(
+    X,
+    parameters,
+    normalise=False,
+    n_jobs=1,
+    original_implementation=False,
+):
+    """Apply a fitted MultiRocket transform.
+
+    Parameters
+    ----------
+    X : 3D np.ndarray, shape (n_instances, n_columns, n_timepoints)
+        panel of time series, made contiguous; float32 is kept and computed in
+        float32, other dtypes are cast to float64. sktime casts all dtypes to
+        float64, so float32 results differ slightly from sktime's.
+    parameters : tuple
+        the tuple returned by ``fit``.
+    normalise : bool, default=False
+        whether to normalise each series before transforming; must match the
+        value used in ``fit``.
+    n_jobs : int, default=1
+        threads for the GIL-releasing Cython kernel over disjoint instance
+        chunks. ``-1`` (or out of range) uses all processors.
+    original_implementation : bool, default=False
+        whether to reproduce sktime's reference behaviour for the differenced
+        pass rather than the corrected transform; see sktime#11291.
+
+    Returns
+    -------
+    np.ndarray, shape (n_instances, n_features), float32
+    """
+    X = np.ascontiguousarray(
+        X, dtype=np.float32 if X.dtype == np.float32 else np.float64
+    )
+    if normalise:
+        X = _normalise(X)
+    X1 = np.ascontiguousarray(np.diff(X, 1))
+
+    _check_dilations(parameters[0][2], X.shape[2], "dilations")
+    _check_dilations(parameters[1][2], X1.shape[2], "dilations1")
+
+    args = (
+        *parameters[0],
+        *parameters[1],
+        original_implementation,
+    )
+    n_instances = X.shape[0]
+
+    if n_jobs < 1 or n_jobs > multiprocessing.cpu_count():
+        n_jobs = multiprocessing.cpu_count()
+    n_jobs = min(n_jobs, n_instances)
+
+    # nan_to_num in place: the kernel's output is freshly allocated, and per
+    # chunk it runs in the worker threads and avoids a second full copy.
+    if n_jobs <= 1:
+        return np.nan_to_num(_cy.transform(X, X1, *args), copy=False)
+
+    # the Cython kernel releases the GIL, so plain threads run truly in
+    # parallel across disjoint instance chunks.
+    bounds = np.linspace(0, n_instances, n_jobs + 1).astype(int)
+    chunks = [
+        (
+            np.ascontiguousarray(X[bounds[i] : bounds[i + 1]]),
+            np.ascontiguousarray(X1[bounds[i] : bounds[i + 1]]),
+        )
+        for i in range(n_jobs)
+        if bounds[i + 1] > bounds[i]
+    ]
+    with ThreadPoolExecutor(max_workers=n_jobs) as ex:
+        parts = list(
+            ex.map(
+                lambda c: np.nan_to_num(_cy.transform(*c, *args), copy=False), chunks
+            )
+        )
+    return np.vstack(parts)
diff --git a/sktime_cython/transformations/rocket/_multirocket_multivariate_cython.pyi b/sktime_cython/transformations/rocket/_multirocket_multivariate_cython.pyi
new file mode 100644
index 0000000..5d622d7
--- /dev/null
+++ b/sktime_cython/transformations/rocket/_multirocket_multivariate_cython.pyi
@@ -0,0 +1,31 @@
+"""Type stubs for the compiled MultiRocketMultivariate Cython kernels."""
+
+import numpy as np
+from numpy.typing import NDArray
+
+_Real = NDArray[np.float32] | NDArray[np.float64]
+
+def transform(
+    X: _Real,
+    X1: _Real,
+    num_channels_per_combination: NDArray[np.int32],
+    channel_indices: NDArray[np.int32],
+    dilations: NDArray[np.int32],
+    num_features_per_dilation: NDArray[np.int32],
+    biases: NDArray[np.float32],
+    num_channels_per_combination1: NDArray[np.int32],
+    channel_indices1: NDArray[np.int32],
+    dilations1: NDArray[np.int32],
+    num_features_per_dilation1: NDArray[np.int32],
+    biases1: NDArray[np.float32],
+    num_features_per_kernel: int,
+    original_implementation: bool,
+) -> NDArray[np.float32]: ...
+def fit_biases(
+    X: _Real,
+    num_channels_per_combination: NDArray[np.int32],
+    channel_indices: NDArray[np.int32],
+    dilations: NDArray[np.int32],
+    num_features_per_dilation: NDArray[np.int32],
+    instance_indices: NDArray[np.int32],
+) -> NDArray[np.float32]: ...
diff --git a/sktime_cython/transformations/rocket/_multirocket_multivariate_cython.pyx b/sktime_cython/transformations/rocket/_multirocket_multivariate_cython.pyx
new file mode 100644
index 0000000..3afb738
--- /dev/null
+++ b/sktime_cython/transformations/rocket/_multirocket_multivariate_cython.pyx
@@ -0,0 +1,431 @@
+# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
+"""Cython MultiRocketMultivariate kernels.
+
+Ahead-of-time compiled ports of sktime's ``_multirocket_multi_numba`` kernels
+(``_fit_biases`` and ``_transform``). Same math, no numba JIT
+warmup. TODO.
+"""
+
+import numpy as np
+
+cimport numpy as cnp
+cimport cython
+from libc.stdlib cimport free, malloc
+from libc.stdint cimport int64_t
+from cython cimport floating
+
+cnp.import_array()
+
+# combinations(range(9), 3) -> 84 kernels, flattened to 252 int32.
+cdef int[252] _IDX
+cdef int _fill_idx():
+    cdef int n = 0, i, j, k
+    for i in range(9):
+        for j in range(i + 1, 9):
+            for k in range(j + 1, 9):
+                _IDX[3 * n] = i
+                _IDX[3 * n + 1] = j
+                _IDX[3 * n + 2] = k
+                n += 1
+    return n
+cdef int _NUM_KERNELS = _fill_idx()
+
+
+cdef void _one(
+    floating[:, :, ::1] X,
+    int ex,
+    int n_columns,
+    int n_timepoints,
+    int raw_n_timepoints,
+    int[::1] num_channels_per_combination,
+    int[::1] channel_indices,
+    int[::1] dilations,
+    int[::1] num_features_per_dilation,
+    float[::1] biases,
+    int num_features,
+    int64_t feature_offset,
+    bint advance_combination,
+    float* C_alpha,
+    float* C_gamma,
+    float* C,
+    float[:, ::1] feat,
+) noexcept nogil:
+    """One representation (raw or differenced) for one instance.
+
+    Writes 4 features per kernel at ``feature_offset + fis + fc`` and every
+    ``num_features`` thereafter. Buffers are caller-allocated and reused.
+    """
+    cdef int num_kernels = _NUM_KERNELS
+    cdef int num_dilations = dilations.shape[0]
+    cdef int csize = n_columns * n_timepoints
+
+    cdef int di, ki, g, c, t, ch
+    cdef int padding, padding0, padding1, nfd, dilation
+    cdef int start, end, base, combination_index
+    cdef int fis, fie, ncc, ncs, nce
+    cdef int i0, i1, i2, fc, n_valid, j
+    cdef float* _c
+    cdef int ppv, last_val, stretch, max_stretch
+    cdef int64_t mean_index, feature_index
+    cdef float bias
+    cdef floating x
+    cdef double mean
+
+    fis = 0
+    ncs = 0
+    combination_index = 0
+    for di in range(num_dilations):
+        padding0 = di % 2
+        dilation = dilations[di]
+        padding = ((9 - 1) * dilation) // 2
+        nfd = num_features_per_dilation[di]
+
+        # C_alpha = A = -X ; C_gamma = 0 ; C_gamma[4] = G = 3X
+        for c in range(n_columns):
+            base = c * n_timepoints
+            for t in range(n_timepoints):
+                x = X[ex, c, t]
+                C_alpha[base + t] = -x
+                C_gamma[4 * csize + base + t] = x + x + x
+        for g in range(9):
+            if g == 4:
+                continue
+            for c in range(n_columns):
+                base = g * csize + c * n_timepoints
+                for t in range(n_timepoints):
+                    C_gamma[base + t] = 0.0
+
+        # gamma_index 0..3. The caller passes raw_n_timepoints == n_timepoints
+        # for the corrected transform; in the original implementation the
+        # differenced pass receives the undifferenced length here, which sits
+        # its windows one sample further left (sktime#11291).
+        end = raw_n_timepoints - padding
+        # `end` stays within the row (and `n_valid` below stays positive) because
+        # dilation >= 1 and 8 * dilation < n_timepoints; boundscheck is off here,
+        # so the Python layer enforces both in ``_common._check_dilations``.
+        for g in range(4):
+            if end > 0:
+                for c in range(n_columns):
+                    base = c * n_timepoints
+                    for t in range(end):
+                        C_alpha[base + n_timepoints - end + t] += -X[ex, c, t]
+                        C_gamma[g * csize + base + n_timepoints - end + t] = (
+                            3.0 * X[ex, c, t]
+                        )
+            end += dilation
+
+        # gamma_index 5..8
+        start = dilation
+        for g in range(5, 9):
+            if start < n_timepoints:
+                for c in range(n_columns):
+                    base = c * n_timepoints
+                    for t in range(n_timepoints - start):
+                        C_alpha[base + t] += -X[ex, c, start + t]
+                        C_gamma[g * csize + base + t] = (
+                            3.0 * X[ex, c, start + t]
+                        )
+            start += dilation
+
+        for ki in range(num_kernels):
+            fie = fis + nfd
+            ncc = num_channels_per_combination[combination_index]
+            nce = ncs + ncc
+            padding1 = (padding0 + ki) % 2
+            i0 = _IDX[3 * ki]
+            i1 = _IDX[3 * ki + 1]
+            i2 = _IDX[3 * ki + 2]
+
+            # C[t] = sum over the combination's channels
+            for t in range(n_timepoints):
+                C[t] = 0.0
+            for ch in range(ncs, nce):
+                c = channel_indices[ch]
+                base = c * n_timepoints
+                for t in range(n_timepoints):
+                    C[t] += (
+                        C_alpha[base + t]
+                        + C_gamma[i0 * csize + base + t]
+                        + C_gamma[i1 * csize + base + t]
+                        + C_gamma[i2 * csize + base + t]
+                    )
+
+            if padding1 == 0:
+                for fc in range(nfd):
+                    bias = biases[fis + fc]
+                    ppv = 0
+                    last_val = 0
+                    max_stretch = 0
+                    mean_index = 0
+                    mean = 0.0
+
+                    for j in range(n_timepoints):
+                        if C[j] > bias:
+                            ppv += 1
+                            mean_index += j
+                            mean += C[j] + bias
+                        elif C[j] < bias:
+                            stretch = j - last_val
+                            if stretch > max_stretch:
+                                max_stretch = stretch
+                            last_val = j
+                    stretch = n_timepoints - 1 - last_val
+                    if stretch > max_stretch:
+                        max_stretch = stretch
+
+                    feature_index = feature_offset + fis + fc
+                    feat[ex, feature_index] = <float>ppv / n_timepoints
+                    feature_index += num_features
+                    feat[ex, feature_index] = <float>max_stretch
+                    feature_index += num_features
+                    feat[ex, feature_index] = mean / ppv if ppv > 0 else 0.0
+                    feature_index += num_features
+                    feat[ex, feature_index] = <float>mean_index / ppv if ppv > 0 else -1.0
+            else:
+                _c = C + padding
+                n_valid = n_timepoints - 2 * padding
+
+                for fc in range(nfd):
+                    bias = biases[fis + fc]
+                    ppv = 0
+                    last_val = 0
+                    max_stretch = 0
+                    mean_index = 0
+                    mean = 0
+
+                    for j in range(n_valid):
+                        if _c[j] > bias:
+                            ppv += 1
+                            mean_index += j
+                            mean += _c[j] + bias
+                        elif _c[j] < bias:
+                            stretch = j - last_val
+                            if stretch > max_stretch:
+                                max_stretch = stretch
+                            last_val = j
+                    stretch = n_valid - 1 - last_val
+                    if stretch > max_stretch:
+                        max_stretch = stretch
+
+                    feature_index = feature_offset + fis + fc
+                    feat[ex, feature_index] = <float>ppv / n_valid
+                    feature_index += num_features
+                    feat[ex, feature_index] = <float>max_stretch
+                    feature_index += num_features
+                    feat[ex, feature_index] = mean / ppv if ppv > 0 else 0.0
+                    feature_index += num_features
+                    feat[ex, feature_index] = (
+                        <float>mean_index / ppv if ppv > 0 else -1.0
+                    )
+
+            fis = fie
+            # the original implementation does not advance these in the
+            # differenced pass, so every kernel there reuses the first
+            # combination's channel selection (sktime#11291).
+            if advance_combination:
+                combination_index += 1
+                ncs = nce
+
+
+def transform(
+    floating[:,:,::1] X,
+    floating[:,:,::1] X1,
+    int[::1] num_channels_per_combination,
+    int[::1] channel_indices,
+    int[::1] dilations,
+    int[::1] num_features_per_dilation,
+    float[::1] biases,
+    int[::1] num_channels_per_combination1,
+    int[::1] channel_indices1,
+    int[::1] dilations1,
+    int[::1] num_features_per_dilation1,
+    float[::1] biases1,
+    bint original_implementation,
+):
+    """Port of ``_transform``. X is (n_instances, n_columns, n_timepoints);
+        X1 is (n_instances, n_columns, n_timepoints-1).
+
+    ``original_implementation`` reproduces sktime's reference behaviour for the
+    differenced pass: windows sized from the undifferenced length, cursors not
+    advanced, and the base pass's channel selection reused.
+    """
+    cdef int n_instances = X.shape[0]
+    cdef int n_columns = X.shape[1]
+    cdef int n_timepoints = X.shape[2]
+
+    cdef int num_kernels = _NUM_KERNELS
+    cdef int num_dilations = dilations.shape[0]
+    cdef int num_dilations1 = dilations1.shape[0]
+
+    cdef int total_fpd = 0
+    cdef int d
+    for d in range(num_dilations):
+        total_fpd += num_features_per_dilation[d]
+    cdef int total_fpd1 = 0
+    cdef int d1
+    for d1 in range(num_dilations1):
+        total_fpd1 += num_features_per_dilation1[d1]
+    cdef int num_features = num_kernels * total_fpd
+    cdef int num_features1 = num_kernels * total_fpd1
+
+    cdef int64_t num_total_features = <int64_t>(num_features + num_features1) * 4
+    cdef cnp.ndarray[cnp.float32_t, ndim=2, mode="c"] features = np.zeros(
+        (n_instances, num_total_features), dtype=np.float32
+    )
+    cdef int64_t num_features_per_transform = num_total_features >> 1
+    cdef float[:, ::1] feat = features
+
+    cdef int csize = n_columns * n_timepoints
+    # Per-instance work buffers (single-threaded; reused across dilations).
+    cdef float* C_alpha = <float*>malloc(csize * sizeof(float))
+    cdef float* C_gamma = <float*>malloc(9 * csize * sizeof(float))
+    cdef float* C = <float*>malloc(n_timepoints * sizeof(float))
+    if C_alpha == NULL or C_gamma == NULL or C == NULL:
+        free(C_alpha); free(C_gamma); free(C)
+        raise MemoryError()
+
+    cdef int ex
+
+    try:
+      with nogil:
+        for ex in range(n_instances):
+            _one(
+                X, ex, n_columns, n_timepoints, n_timepoints,
+                num_channels_per_combination, channel_indices,
+                dilations, num_features_per_dilation, biases,
+                num_features, 0, True,
+                C_alpha, C_gamma, C, feat,
+            )
+            if original_implementation:
+                _one(
+                    X1, ex, n_columns, n_timepoints - 1, n_timepoints,
+                    num_channels_per_combination, channel_indices,
+                    dilations1, num_features_per_dilation1, biases1,
+                    num_features, num_features_per_transform, False,
+                    C_alpha, C_gamma, C, feat,
+                )
+            else:
+                _one(
+                    X1, ex, n_columns, n_timepoints - 1, n_timepoints - 1,
+                    num_channels_per_combination1, channel_indices1,
+                    dilations1, num_features_per_dilation1, biases1,
+                    num_features, num_features_per_transform, True,
+                    C_alpha, C_gamma, C, feat,
+                )
+
+    finally:
+        free(C_alpha)
+        free(C_gamma)
+        free(C)
+
+    return features
+
+
+def fit_biases(
+    floating[:,:,::1] X,
+    int[::1] num_channels_per_combination,
+    int[::1] channel_indices,
+    int[::1] dilations,
+    int[::1] num_features_per_dilation,
+    int[::1] instance_indices,
+):
+    """Build per-combination convolution output C, summed over channels.
+
+    Port of the convolution-building half of ``_fit_biases``. Returns a
+    ``(num_combinations, n_timepoints)`` float32 array; the caller applies
+    ``np.quantile`` per combination to obtain biases. ``instance_indices`` holds
+    the ``np.random.randint(n_instances)`` draw for each combination, computed by
+    the caller to reproduce the numba random sequence exactly.
+    """
+    cdef int n_timepoints = X.shape[2]
+
+    cdef int num_kernels = _NUM_KERNELS
+    cdef int num_dilations = dilations.shape[0]
+    cdef int num_combinations = num_kernels * num_dilations
+
+    cdef cnp.ndarray[cnp.float32_t, ndim=2, mode="c"] out = np.zeros(
+        (num_combinations, n_timepoints), dtype=np.float32
+    )
+    cdef float[:, ::1] outv = out
+
+    # 9 gamma planes x n_timepoints for a single (instance, channel) row.
+    cdef float* C_alpha = <float*>malloc(n_timepoints * sizeof(float))
+    cdef float* C_gamma = <float*>malloc(9 * n_timepoints * sizeof(float))
+    if C_alpha == NULL or C_gamma == NULL:
+        free(C_alpha); free(C_gamma)
+        raise MemoryError()
+
+    cdef int di, ki, g, t, ch, ex, c
+    cdef int padding, dilation, start, end
+    cdef int comb, ncs, nce, ncc
+    cdef int i0, i1, i2
+    cdef floating x
+
+    try:
+        comb = 0
+        ncs = 0
+        for di in range(num_dilations):
+            dilation = dilations[di]
+            padding = ((9 - 1) * dilation) // 2
+
+            for ki in range(num_kernels):
+                ncc = num_channels_per_combination[comb]
+                nce = ncs + ncc
+                ex = instance_indices[comb]
+                i0 = _IDX[3 * ki]
+                i1 = _IDX[3 * ki + 1]
+                i2 = _IDX[3 * ki + 2]
+
+                for t in range(n_timepoints):
+                    outv[comb, t] = 0.0
+
+                # Accumulate C = C_alpha + C_gamma[i0,i1,i2], summed over channels.
+                for ch in range(ncs, nce):
+                    c = channel_indices[ch]
+
+                    # C_alpha = A = -X ; C_gamma = 0 ; C_gamma[4] = G = 3X
+                    for t in range(n_timepoints):
+                        x = X[ex, c, t]
+                        C_alpha[t] = -x
+                        C_gamma[4 * n_timepoints + t] = x + x + x
+                    for g in range(9):
+                        if g == 4:
+                            continue
+                        for t in range(n_timepoints):
+                            C_gamma[g * n_timepoints + t] = 0.0
+
+                    end = n_timepoints - padding
+                    for g in range(4):
+                        if end > 0:
+                            for t in range(end):
+                                C_alpha[n_timepoints - end + t] += -X[ex, c, t]
+                                C_gamma[g * n_timepoints + n_timepoints - end + t] = (
+                                    3.0 * X[ex, c, t]
+                                )
+                        end += dilation
+
+                    start = dilation
+                    for g in range(5, 9):
+                        if start < n_timepoints:
+                            for t in range(n_timepoints - start):
+                                C_alpha[t] += -X[ex, c, start + t]
+                                C_gamma[g * n_timepoints + t] = (
+                                    3.0 * X[ex, c, start + t]
+                                )
+                        start += dilation
+
+                    for t in range(n_timepoints):
+                        outv[comb, t] += (
+                            C_alpha[t]
+                            + C_gamma[i0 * n_timepoints + t]
+                            + C_gamma[i1 * n_timepoints + t]
+                            + C_gamma[i2 * n_timepoints + t]
+                        )
+
+                comb += 1
+                ncs = nce
+    finally:
+        free(C_alpha)
+        free(C_gamma)
+
+    return out
diff --git a/sktime_cython/transformations/rocket/tests/test_common.py b/sktime_cython/transformations/rocket/tests/test_common.py
new file mode 100644
index 0000000..f60a32e
--- /dev/null
+++ b/sktime_cython/transformations/rocket/tests/test_common.py
@@ -0,0 +1,72 @@
+"""Tests for the scaffolding shared by MiniRocket and MultiRocket.
+
+Self-contained (numpy + pytest only) so they run in cibuildwheel's isolated
+wheel-test env.
+"""
+
+import numpy as np
+import pytest
+
+from sktime_cython.transformations.rocket._common import (
+    _biases_from_C,
+    _fit_dilations,
+    _quantiles,
+)
+
+
+def _biases_from_C_loop(C, quantiles, num_features_per_dilation, num_kernels):
+    """Reference: numba's np.quantile, per combination, in a plain loop.
+
+    numba/np/arraymath.py promotes to float64, scales the quantile to a
+    percentile and back, and interpolates between ranks f-1 and f. Spelled out
+    here so the vectorized version can be checked against something readable,
+    without needing numba installed.
+    """
+    biases = np.zeros(
+        num_kernels * int(np.sum(num_features_per_dilation)), dtype=np.float32
+    )
+    feature_index_start = 0
+    combination_index = 0
+    for dilation_index in range(len(num_features_per_dilation)):
+        nfd = num_features_per_dilation[dilation_index]
+        for _kernel_index in range(num_kernels):
+            feature_index_end = feature_index_start + nfd
+            a = np.sort(C[combination_index]).astype(np.float64)
+            for i in range(feature_index_start, feature_index_end):
+                percentile = np.float64(quantiles[i]) * 100.0
+                rank = 1.0 + (len(a) - 1) * (percentile / 100.0)
+                floor = np.floor(rank)
+                weight = rank - floor
+                lower = int(floor) - 1
+                biases[i] = np.float32(
+                    a[lower] * (1 - weight) + a[min(lower + 1, len(a) - 1)] * weight
+                )
+            feature_index_start = feature_index_end
+            combination_index += 1
+    return biases
+
+
+@pytest.mark.parametrize("n_timepoints", [9, 60, 137])
+@pytest.mark.parametrize(
+    "num_features_per_dilation", [[1], [4, 2, 1], [32, 16, 8, 4, 2, 1, 1]]
+)
+def test_biases_from_C_matches_loop(num_features_per_dilation, n_timepoints):
+    """Vectorized quantiles must match the per-combination loop exactly."""
+    num_kernels = 84
+    nfpd = np.array(num_features_per_dilation, dtype=np.int32)
+    rng = np.random.RandomState(0)
+    C = rng.normal(size=(num_kernels * len(nfpd), n_timepoints)).astype(np.float32)
+    quantiles = _quantiles(num_kernels * int(nfpd.sum()))
+
+    np.testing.assert_array_equal(
+        _biases_from_C(C, quantiles, nfpd, num_kernels),
+        _biases_from_C_loop(C, quantiles, nfpd, num_kernels),
+    )
+
+
+@pytest.mark.parametrize("num_features", [84, 168, 10_000])
+def test_fit_dilations_features_add_up(num_features):
+    """Features per dilation sum to the per-kernel budget, one per dilation."""
+    dilations, nfpd = _fit_dilations(137, num_features, 32)
+    assert len(dilations) == len(nfpd)
+    assert nfpd.sum() == max(num_features, 84) // 84
diff --git a/sktime_cython/transformations/rocket/tests/test_multirocket_cython.py b/sktime_cython/transformations/rocket/tests/test_multirocket_cython.py
new file mode 100644
index 0000000..01e520e
--- /dev/null
+++ b/sktime_cython/transformations/rocket/tests/test_multirocket_cython.py
@@ -0,0 +1,132 @@
+"""Tests for the Cython MultiRocket transform.
+
+Mirrors ``test_minirocket_cython``: the shape / threading / guard tests are
+self-contained (no sktime) so they run in cibuildwheel's isolated wheel-test
+env, and the equivalence-vs-numba test is skipped where sktime is absent.
+"""
+
+import importlib.util
+
+import numpy as np
+import pytest
+
+from sktime_cython.transformations.rocket._multirocket import (
+    multirocket_fit,
+    multirocket_transform,
+)
+
+_HAS_SKTIME = importlib.util.find_spec("sktime") is not None
+
+
+def _sktime_has_flag():
+    """Whether the installed sktime exposes ``original_implementation``."""
+    if not _HAS_SKTIME:
+        return False
+    from sktime.transformations.rocket import MultiRocketMultivariate
+
+    return "original_implementation" in MultiRocketMultivariate.get_param_names()
+
+
+# At large num_kernels a value can sit on the bias threshold and be classified
+# differently by the two float32 summation orders, flipping the PPV/MPV/MIPV
+# triplet of one feature. Measured rate is ~2e-5 of elements; allow 1e-4.
+_MISMATCH_BUDGET = 1e-4
+
+
+def _panel(seed, n_columns=3, n_timepoints=60):
+    rng = np.random.RandomState(seed)
+    return rng.normal(size=(6, n_columns, n_timepoints)).astype(np.float64)
+
+
+@pytest.mark.skipif(not _HAS_SKTIME, reason="sktime not installed (dev extra)")
+@pytest.mark.parametrize("original_implementation", [False, True])
+@pytest.mark.parametrize("n_columns", [1, 4])
+@pytest.mark.parametrize(
+    "num_kernels,max_dilations_per_kernel,random_state,n_timepoints",
+    [(84, 32, 42, 60), (168, 16, 7, 60), (6250, 32, 0, 137)],
+)
+def test_cython_matches_numba(
+    n_columns,
+    num_kernels,
+    max_dilations_per_kernel,
+    random_state,
+    n_timepoints,
+    original_implementation,
+):
+    """Cython transform must match the numba implementation (groundtruth)."""
+    from sktime.transformations.rocket import MultiRocketMultivariate
+
+    kwargs = {}
+    if _sktime_has_flag():
+        kwargs["original_implementation"] = original_implementation
+    elif not original_implementation:
+        # sktime before the fix only implements the original behaviour
+        pytest.skip("installed sktime predates original_implementation")
+
+    X = _panel(random_state, n_columns=n_columns, n_timepoints=n_timepoints)
+
+    numba_out = MultiRocketMultivariate(
+        num_kernels=num_kernels,
+        max_dilations_per_kernel=max_dilations_per_kernel,
+        random_state=random_state,
+        **kwargs,
+    ).fit_transform(X)
+
+    params = multirocket_fit(
+        X,
+        num_kernels=num_kernels,
+        max_dilations_per_kernel=max_dilations_per_kernel,
+        random_state=random_state,
+    )
+    cython_out = multirocket_transform(
+        X, params, original_implementation=original_implementation
+    )
+
+    assert cython_out.shape == numba_out.shape
+    close = np.isclose(cython_out, numba_out.to_numpy(), rtol=1e-4, atol=1e-5)
+    assert 1 - close.mean() <= _MISMATCH_BUDGET
+
+
+@pytest.mark.parametrize("num_kernels", [84, 168])
+def test_output_shape(num_kernels):
+    """transform yields 4 features per kernel, on X and its difference."""
+    X = _panel(0)
+    out = multirocket_transform(
+        X, multirocket_fit(X, num_kernels=num_kernels, random_state=0)
+    )
+    assert out.shape == (X.shape[0], 2 * 4 * (num_kernels // 84) * 84)
+    assert out.dtype == np.float32
+
+
+def test_threaded_matches_serial():
+    """n_jobs>1 must match the single-threaded result."""
+    X = _panel(3)
+    params = multirocket_fit(X, num_kernels=168, random_state=1)
+    np.testing.assert_array_equal(
+        multirocket_transform(X, params, n_jobs=1),
+        multirocket_transform(X, params, n_jobs=4),
+    )
+
+
+def test_float32_and_float64_both_accepted():
+    """The kernels are fused over float32/float64; both dispatch."""
+    X = _panel(5)
+    params = multirocket_fit(X, num_kernels=84, random_state=0)
+    out64 = multirocket_transform(X, params)
+    out32 = multirocket_transform(X.astype(np.float32), params)
+    assert out64.shape == out32.shape
+    # the two dtypes round differently, so a value sitting on a bias threshold
+    # can be classified either way; budget for that rather than requiring
+    # elementwise closeness
+    close = np.isclose(out64, out32, rtol=1e-3, atol=1e-3)
+    assert 1 - close.mean() <= 1e-2
+
+
+@pytest.mark.parametrize("bad", [0, 40])
+def test_invalid_dilations_raise(bad):
+    """Dilations outside the kernel's in-bounds invariants are rejected."""
+    X = _panel(0)
+    p, p1 = multirocket_fit(X, num_kernels=84, random_state=0)
+    p = (*p[:2], np.array([bad], dtype=np.int32), *p[3:])
+    with pytest.raises(ValueError, match="dilations"):
+        multirocket_transform(X, (p, p1))