diff --git a/.github/workflows/test_and_deploy.yml b/.github/workflows/test_and_deploy.yml
index 8916816..befa33e 100644
--- a/.github/workflows/test_and_deploy.yml
+++ b/.github/workflows/test_and_deploy.yml
@@ -24,7 +24,7 @@ jobs:
     strategy:
       matrix:
         platform: [ubuntu-latest, windows-latest, macos-latest]
-        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13']
+        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13', '3.14']
 
     steps:
       - uses: actions/checkout@v4
@@ -43,11 +43,18 @@ jobs:
           python -m pip install --upgrade pip
           python -m pip install setuptools tox tox-gh-actions
 
+      - name: Set up MSVC
+        if: runner.os == 'Windows'
+        uses: TheMrMilchmann/setup-msvc-dev@v4
+        with:
+          arch: x64
+
       # this runs the platform-specific tests declared in tox.ini
       - name: Test with tox
-        run: python -m tox
         env:
           PLATFORM: ${{ matrix.platform }}
+          TORCHINDUCTOR_CACHE_DIR: ${{ runner.temp }}/torchinductor
+        run: python -m tox
 
       - name: Upload coverage reports to Codecov
         uses: codecov/codecov-action@v4
diff --git a/examples/skip_performance_basis_functions.py b/examples/skip_performance_basis_functions.py
new file mode 100644
index 0000000..f43f058
--- /dev/null
+++ b/examples/skip_performance_basis_functions.py
@@ -0,0 +1,54 @@
+import collections
+import time
+
+import jax.numpy as jax
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import seaborn as sns
+import splinebox
+import torch
+
+results = collections.defaultdict(list)
+
+for length in [10, 100, 1000, 10000, 100000]:
+    for replicate in np.arange(11):
+        x = np.random.rand(length)
+        for array_namespace in ["numpy", "torch", "jax"]:
+            if array_namespace == "torch":
+                x = torch.tensor(x)
+            elif array_namespace == "jax":
+                x = jax.asarray(x)
+            for basis_function_name, basis_function_class in splinebox.basis_functions.inventory().items():
+                if "Exponential" in basis_function_name:
+                    basis_function = basis_function_class(5)
+                else:
+                    basis_function = basis_function_class()
+                for deriv in range(3):
+                    for jit in [True, False]:
+                        start = time.perf_counter_ns()
+                        basis_function(x, derivative=deriv, jit=jit)
+                        stop = time.perf_counter_ns()
+
+                        results["length"].append(length)
+                        results["Replicate"].append(replicate)
+                        # results["Array namespace"].append(array_namespace)
+                        results["Basis function"].append(basis_function_name)
+                        results["Derivative"].append(deriv)
+                        results["JIT"].append(jit)
+                        results["time [ns]"].append(stop - start)
+
+df = pd.DataFrame(results)
+
+g = sns.relplot(
+    data=df.loc[df["Replicate"] > 0],
+    x="length",
+    y="time [ns]",
+    col="Basis function",
+    hue="JIT",
+    row="Array namespace",
+    kind="line",
+)
+g.set(xscale="log")
+g.set(yscale="log")
+plt.show()
diff --git a/pyproject.toml b/pyproject.toml
index 794f5d8..c6fa985 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -94,6 +94,9 @@ test = [
     "tox",
     "pytest",
     "pytest-cov",
+    "array-api-compat",
+    "torch",
+    "jax"
 ]
 examples = [
     "matplotlib",
diff --git a/src/splinebox/basis_functions.py b/src/splinebox/basis_functions.py
index d69a590..fa2d8b1 100644
--- a/src/splinebox/basis_functions.py
+++ b/src/splinebox/basis_functions.py
@@ -10,15 +10,107 @@
 For more information on implementing a new basis function, see :class:`splinebox.basis_functions.BasisFunction`.
 """
 
+import collections
+import contextlib
 import inspect
 import math
+import numbers
 import sys
 import warnings
 
+import array_api_compat
+import array_api_compat.numpy
 import numba
 import numpy as np
 
 
+@contextlib.contextmanager
+def _torch_compile_recompile_limit(limit=64):
+    """
+    Context manager that temporarily raises the torch dynamo recompile limit.
+
+    Compiling the filters with torch.compile unrolls their loops statically,
+    which means a separate graph is compiled for every combination of input
+    shape and basis function parameters (e.g. the pole of the exponential
+    basis function depends on ``M``). This can exceed the default recompile
+    limit of a single code object, so it is temporarily raised while compiling.
+    """
+    import torch
+
+    config = torch._dynamo.config
+    key = "recompile_limit" if hasattr(config, "recompile_limit") else "cache_size_limit"
+    old = getattr(config, key)
+    setattr(config, key, max(old, limit))
+    try:
+        yield
+    finally:
+        setattr(config, key, old)
+
+
+def _cached_compiled(fn, key, compile_fn):
+    """
+    Return a compiled version of ``fn``, caching the compiled wrapper.
+
+    Constructing compiled wrappers, e.g. with ``torch.compile`` or
+    ``jax.jit``, is expensive, so each wrapper is constructed only once.
+    The wrappers are cached in the ``_compiled_wrapper_cache`` dictionary
+    of the owning class. Wrappers compiled for bound methods capture the
+    instance they were compiled for, so the cache is keyed by the
+    instance's attributes in addition to the function. This allows fresh
+    instances with the same state, e.g. a new ``Exponential(5)``, to
+    reuse the wrapper compiled for a previous one, while instances with
+    different state, e.g. ``Exponential(7)``, get their own wrapper. It
+    relies on basis function attributes not being modified after the
+    first compiled call, which holds because they are only set in
+    ``__init__``. Note that this keeps the first instance with a given
+    state alive for the lifetime of the process. Wrappers for plain
+    functions are keyed by the function itself. ``key`` identifies the
+    compilation configuration and must be distinct for functions that are
+    compiled differently, e.g. for different backends.
+    """
+    if hasattr(fn, "__self__"):
+        instance = fn.__self__
+        owner = type(instance)
+        cache_key = (fn.__func__, key, tuple(sorted(instance.__dict__.items())))
+    else:
+        parts = fn.__qualname__.split(".")[:-1]  # drop the function name
+        owner = sys.modules[fn.__module__]
+        for part in parts:
+            owner = getattr(owner, part)
+        cache_key = (fn, key)
+    if cache_key not in owner._compiled_wrapper_cache:
+        owner._compiled_wrapper_cache[cache_key] = compile_fn(fn)
+    return owner._compiled_wrapper_cache[cache_key]
+
+
+def _item_set(arr, idx, value, xp):
+    """
+    Assign ``value`` to ``arr[idx]`` and return the updated array.
+
+    JAX arrays are immutable and do not support in-place item assignment,
+    so the assignment is performed functionally using ``arr.at[idx].set``.
+    All other backends support in-place assignment directly.
+    """
+    if array_api_compat.is_jax_namespace(xp):
+        return arr.at[idx].set(value)
+    arr[idx] = value
+    return arr
+
+
+def _item_add(arr, idx, value, xp):
+    """
+    Add ``value`` to ``arr[idx]`` and return the updated array.
+
+    JAX arrays are immutable and do not support in-place item assignment,
+    so the addition is performed functionally using ``arr.at[idx].add``.
+    All other backends support in-place assignment directly.
+    """
+    if array_api_compat.is_jax_namespace(xp):
+        return arr.at[idx].add(value)
+    arr[idx] += value
+    return arr
+
+
 class BasisFunction:
     """
     Base class for all basis functions.
@@ -54,15 +146,15 @@ class BasisFunction:
     ...         # Change this if your new basis function is not in splinebox.basis_functions
     ...         return "splinebox.basis_functions.MyBasis()"
     ...
-    ...     def _func(self, t):
+    ...     def _func_numba(self, t):
     ...         # Implement your function here
     ...         return val
     ...
-    ...     def _derivative_1(self, t):
+    ...     def _derivative_1_numba(self, t):
     ...         # Implement the first derivative of your function here
     ...         return val
     ...
-    ...     def _derivative_2(self, t):
+    ...     def _derivative_2_numba(self, t):
     ...         # Implement the second derivative of your function here
     ...         # or raise an error if not differentiable
     ...         raise RuntimeError("MyBasis isn't twice differentiable.")
@@ -79,6 +171,7 @@ class BasisFunction:
     """
 
     _unimplemented_message = "This function is not implemented."
+    _compiled_wrapper_cache: dict[tuple, collections.abc.Callable] = {}
 
     def __init__(self, multigenerator, support):
         self.multigenerator = multigenerator
@@ -97,7 +190,7 @@ def __eq__(self, other):
             and other.support == self.support
         )
 
-    def __call__(self, t, derivative=0):
+    def __call__(self, t, derivative=0, jit=True):
         """
         Evaluate the function at position(s) t. You can optionally
         pass derivative=1 or 2 to obtain the respective derivative(s).
@@ -128,25 +221,55 @@ def __call__(self, t, derivative=0):
         >>> basis_function(t, derivative=1)
         array([ 1., -1.])
         """
+        xp = array_api_compat.numpy if isinstance(t, (int, float)) else array_api_compat.array_namespace(t)
+
+        # The xp functions can only accept arrays to for single value we use the numba version
+        use_numba = (jit and array_api_compat.is_numpy_namespace(xp)) or isinstance(t, numbers.Number)
+
+        # Select the function that should be run
         if derivative == 0:
-            return self._func(t)
+            func = self._func_numba if use_numba else self._func_xp
         elif derivative == 1:
-            return self._derivative_1(t)
+            func = self._derivative_1_numba if use_numba else self._derivative_1_xp
         elif derivative == 2:
-            return self._derivative_2(t)
+            func = self._derivative_2_numba if use_numba else self._derivative_2_xp
         else:
             raise ValueError(f"derivative has to be 0, 1, or 2 not {derivative}")
 
-    def _func(self, t):
+        # Decide how to run the function
+        if use_numba:
+            return func(t)
+        elif jit and array_api_compat.is_torch_namespace(xp):
+            with _torch_compile_recompile_limit():
+                return _cached_compiled(
+                    func, "torch", lambda fn: xp.compile(fullgraph=True, dynamic=True, mode="reduce-overhead")(fn)
+                )(t, xp)
+        elif jit and array_api_compat.is_jax_namespace(xp):
+            import jax
+
+            return _cached_compiled(func, "jax", lambda fn: jax.jit(fn, static_argnames=["xp"]))(t, xp)
+        else:
+            return func(t, xp)
+
+    def _func_numba(self, t):
+        raise NotImplementedError(BasisFunction._unimplemented_message)
+
+    def _func_xp(self, t, xp):
+        raise NotImplementedError(BasisFunction._unimplemented_message)
+
+    def _derivative_1_numba(self, t):
+        raise NotImplementedError(BasisFunction._unimplemented_message)
+
+    def _derivative_1_xp(self, t, xp):
         raise NotImplementedError(BasisFunction._unimplemented_message)
 
-    def _derivative_1(self, t):
+    def _derivative_2_numba(self, t):
         raise NotImplementedError(BasisFunction._unimplemented_message)
 
-    def _derivative_2(self, t):
+    def _derivative_2_xp(self, t, xp):
         raise NotImplementedError(BasisFunction._unimplemented_message)
 
-    def filter_symmetric(self, s):
+    def filter_symmetric(self, s, jit=False):
         """
         Returns a filtered version of the input s,
         used to convert knots into control points for an **open spline**.
@@ -191,9 +314,46 @@ def filter_symmetric(self, s):
 
         To perform the same test for the first and the last knot, padding is required.
         """
+        xp = array_api_compat.array_namespace(s)
+
+        if not xp.isdtype(s.dtype, kind="real floating"):
+            # Convert to default floating type
+            s = xp.astype(s, xp.asarray(1.0).dtype)
+
+        use_numba = jit and array_api_compat.is_numpy_namespace(xp)
+
+        if use_numba:
+            result = self._filter_symmetric_numba(s)
+        elif jit and array_api_compat.is_torch_namespace(xp):
+            with _torch_compile_recompile_limit():
+                result = _cached_compiled(
+                    self._filter_symmetric_xp,
+                    "torch",
+                    lambda fn: xp.compile(fullgraph=True, dynamic=False, mode="reduce-overhead")(fn),
+                )(s, xp)
+        elif jit and array_api_compat.is_jax_namespace(xp):
+            import jax
+
+            result = _cached_compiled(self._filter_symmetric_xp, "jax", lambda fn: jax.jit(fn, static_argnames=["xp"]))(
+                s, xp
+            )
+        else:
+            result = self._filter_symmetric_xp(s, xp)
+
+        if s.ndim == 1:
+            # array_api_compat.torch requires the axis argument for squeeze,
+            # so the unit axes are determined explicitly.
+            axes = tuple(axis for axis, size in enumerate(result.shape) if size == 1)
+            return xp.squeeze(result, axis=axes) if axes else result
+        return result
+
+    def _filter_symmetric_numba(self, s):
+        raise NotImplementedError(BasisFunction._unimplemented_message)
+
+    def _filter_symmetric_xp(self, s, xp):
         raise NotImplementedError(BasisFunction._unimplemented_message)
 
-    def filter_periodic(self, s):
+    def filter_periodic(self, s, jit=False):
         """
         Returns a filtered version of the input s, used to convert
         knots into control points for a **closed spline**.
@@ -244,6 +404,43 @@ def filter_periodic(self, s):
         >>> b3(-1) * control_points[1] + b3(0) * control_points[2] + b3(1) * control_points[0]
         array([1., 3.])
         """
+        xp = array_api_compat.array_namespace(s)
+
+        if not xp.isdtype(s.dtype, kind="real floating"):
+            # Convert to default floating type
+            s = xp.astype(s, xp.asarray(1.0).dtype)
+
+        use_numba = jit and array_api_compat.is_numpy_namespace(xp)
+
+        if use_numba:
+            result = self._filter_periodic_numba(s)
+        elif jit and array_api_compat.is_torch_namespace(xp):
+            with _torch_compile_recompile_limit():
+                result = _cached_compiled(
+                    self._filter_periodic_xp,
+                    "torch",
+                    lambda fn: xp.compile(fullgraph=True, dynamic=False, mode="reduce-overhead")(fn),
+                )(s, xp)
+        elif jit and array_api_compat.is_jax_namespace(xp):
+            import jax
+
+            result = _cached_compiled(self._filter_periodic_xp, "jax", lambda fn: jax.jit(fn, static_argnames=["xp"]))(
+                s, xp
+            )
+        else:
+            result = self._filter_periodic_xp(s, xp)
+
+        if s.ndim == 1:
+            # array_api_compat.torch requires the axis argument for squeeze,
+            # so the unit axes are determined explicitly.
+            axes = tuple(axis for axis, size in enumerate(result.shape) if size == 1)
+            return xp.squeeze(result, axis=axes) if axes else result
+        return result
+
+    def _filter_periodic_numba(self, s):
+        raise NotImplementedError(BasisFunction._unimplemented_message)
+
+    def _filter_periodic_xp(self, s, xp):
         raise NotImplementedError(BasisFunction._unimplemented_message)
 
     def refinement_mask(self):
@@ -311,16 +508,24 @@ def __repr__(self):
         return "splinebox.basis_functions.B1()"
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _func(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _func_numba(t):  # pragma: no cover
         val = 0
         if abs(t) >= 0 and abs(t) < 1:
             val = 1 - abs(t)
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _derivative_1(t):  # pragma: no cover
+    def _func_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            xp.abs(t) < 1,
+            1.0 - xp.abs(t),
+            0.0,
+        )
+
+    @staticmethod
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _derivative_1_numba(t):  # pragma: no cover
         val = 0
         if t > -1 and t < 0:
             val = 1
@@ -331,17 +536,36 @@ def _derivative_1(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _derivative_2(t):
+    def _derivative_1_xp(t, xp):  # pragma: no cover
+        abs_t = xp.abs(t)
+        result = -xp.sign(t)
+        result = xp.where((t == 0) | (abs_t == 1), xp.nan, result)
+        return xp.where(abs_t > 1, 0.0, result)
+
+    @staticmethod
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _derivative_2_numba(t):
         return np.nan if t in (-1, 0, 1) else 0
 
     @staticmethod
-    def filter_symmetric(s):
-        return s.astype(float)
+    def _derivative_2_xp(t, xp):
+        return xp.where((t == -1) | (t == 0) | (t == 1), xp.nan, xp.asarray(0.0, dtype=t.dtype))
+
+    @staticmethod
+    def _filter_symmetric_numba(s):
+        return s
+
+    @staticmethod
+    def _filter_symmetric_xp(s, xp):
+        return s
 
     @staticmethod
-    def filter_periodic(s):
-        return s.astype(float)
+    def _filter_periodic_numba(s):
+        return s
+
+    @staticmethod
+    def _filter_periodic_xp(s, xp):
+        return s
 
     def refinement_mask(self):
         order = int(self.support)
@@ -402,8 +626,8 @@ def __repr__(self):
         return "splinebox.basis_functions.B2()"
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _func(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _func_numba(t):  # pragma: no cover
         val = 0
         if t >= -1.5 and t <= -0.5:
             val = 0.5 * (t**2) + 1.5 * t + 1.125
@@ -414,8 +638,8 @@ def _func(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _derivative_1(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _derivative_1_numba(t):  # pragma: no cover
         val = 0
         if t >= -1.5 and t <= -0.5:
             val = t + 1.5
@@ -426,8 +650,8 @@ def _derivative_1(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _derivative_2(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _derivative_2_numba(t):  # pragma: no cover
         val = np.nan
         if t < -1.5:
             val = 0
@@ -441,6 +665,62 @@ def _derivative_2(t):  # pragma: no cover
             val = 0
         return val
 
+    @staticmethod
+    def _func_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= -1.5) & (t <= -0.5),
+            0.5 * (t**2) + 1.5 * t + 1.125,
+            xp.where(
+                (t > -0.5) & (t <= 0.5),
+                -t * t + 0.75,
+                xp.where(
+                    (t > 0.5) & (t <= 1.5),
+                    0.5 * (t**2) - 1.5 * t + 1.125,
+                    0.0,
+                ),
+            ),
+        )
+
+    @staticmethod
+    def _derivative_1_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= -1.5) & (t <= -0.5),
+            t + 1.5,
+            xp.where(
+                (t > -0.5) & (t <= 0.5),
+                -2 * t,
+                xp.where(
+                    (t > 0.5) & (t <= 1.5),
+                    t - 1.5,
+                    0.0,
+                ),
+            ),
+        )
+
+    @staticmethod
+    def _derivative_2_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            t < -1.5,
+            0.0,
+            xp.where(
+                (t > -1.5) & (t < -0.5),
+                1.0,
+                xp.where(
+                    (t > -0.5) & (t < 0.5),
+                    -2.0,
+                    xp.where(
+                        (t > 0.5) & (t < 1.5),
+                        1.0,
+                        xp.where(
+                            t > 1.5,
+                            xp.asarray(0.0, dtype=t.dtype),
+                            xp.nan,
+                        ),
+                    ),
+                ),
+            ),
+        )
+
     def refinement_mask(self):
         order = int(self.support)
         mask = np.zeros(order + 1)
@@ -500,8 +780,8 @@ def __repr__(self):
         return "splinebox.basis_functions.B3()"
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _func(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _func_numba(t):  # pragma: no cover
         val = 0
         if abs(t) >= 0 and abs(t) < 1:
             val = 2 / 3 - (abs(t) ** 2) + (abs(t) ** 3) / 2
@@ -510,8 +790,8 @@ def _func(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _derivative_1(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _derivative_1_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t < 1:
             val = -2 * t + 1.5 * t * t
@@ -524,8 +804,8 @@ def _derivative_1(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _derivative_2(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _derivative_2_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t < 1:
             val = -2 + 3 * t
@@ -537,16 +817,62 @@ def _derivative_2(t):  # pragma: no cover
             val = 2 + t
         return val
 
-    def filter_symmetric(self, s):
-        # This is necessary for the docs to build correctly
-        result = self._filter_symmetric(s.astype(float))
-        if s.ndim == 1:
-            return np.squeeze(result)
-        return result
+    @staticmethod
+    def _func_xp(t, xp):  # pragma: no cover
+        abs_t = xp.abs(t)
+        return xp.where(
+            abs_t < 1,
+            2 / 3 - (abs_t**2) + (abs_t**3) / 2,
+            xp.where(
+                (abs_t >= 1) & (abs_t <= 2),
+                ((2 - abs_t) ** 3) / 6,
+                0.0,
+            ),
+        )
+
+    @staticmethod
+    def _derivative_1_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= 0) & (t < 1),
+            -2 * t + 1.5 * t * t,
+            xp.where(
+                (t > -1) & (t < 0),
+                -2 * t - 1.5 * t * t,
+                xp.where(
+                    (t >= 1) & (t <= 2),
+                    -0.5 * ((2 - t) ** 2),
+                    xp.where(
+                        (t >= -2) & (t <= -1),
+                        0.5 * ((2 + t) ** 2),
+                        0.0,
+                    ),
+                ),
+            ),
+        )
+
+    @staticmethod
+    def _derivative_2_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= 0) & (t < 1),
+            -2 + 3 * t,
+            xp.where(
+                (t > -1) & (t < 0),
+                -2 - 3 * t,
+                xp.where(
+                    (t >= 1) & (t <= 2),
+                    2 - t,
+                    xp.where(
+                        (t >= -2) & (t <= -1),
+                        2 + t,
+                        0.0,
+                    ),
+                ),
+            ),
+        )
 
     @staticmethod
     @numba.jit(nopython=True, nogil=True, cache=True)
-    def _filter_symmetric(s):  # pragma: no cover
+    def _filter_symmetric_numba(s):  # pragma: no cover
         M = len(s)
         pole = -2 + np.sqrt(3)
 
@@ -579,22 +905,48 @@ def _filter_symmetric(s):  # pragma: no cover
         c = c.reshape(shape)
         return c
 
-    def filter_periodic(self, s):
-        # This is necessary for the docs to build correctly
-        result = self._filter_periodic(s.astype(float))
-        if s.ndim == 1:
-            return np.squeeze(result)
-        return result
+    @staticmethod
+    def _filter_symmetric_xp(s, xp):  # pragma: no cover
+        M = len(s)
+        pole = -2 + math.sqrt(3)
+
+        ndim = 1 if s.ndim == 1 else s.shape[1]
+
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(s.dtype, kind="integral") else s.dtype
+        device = s.device if hasattr(s, "device") else None
+        cp = xp.zeros((M, ndim), dtype=dtype, device=device)
+        eps = 1e-8
+        k0 = min(((2 * M) - 2, int(math.ceil(math.log(eps) / math.log(abs(pole))))))
+        for k in range(k0):
+            m = k % (2 * M - 2)
+            val = s[2 * M - 2 - m] if m >= M else s[m]
+            cp = _item_add(cp, 0, val * (pole**m), xp)
+        cp = _item_set(cp, 0, cp[0] * (1 / (1 - (pole ** (2 * M - 2)))), xp)
+
+        for k in range(1, M):
+            cp = _item_set(cp, k, s[k] + pole * cp[k - 1], xp)
+
+        cm = xp.zeros((M, ndim), dtype=dtype, device=device)
+        cm = _item_set(cm, M - 1, cp[M - 1] + (pole * cp[M - 2]), xp)
+        cm = _item_set(cm, M - 1, cm[M - 1] * (pole / ((pole**2) - 1)), xp)
+
+        for k in range(M - 2, -1, -1):
+            cm = _item_set(cm, k, pole * (cm[k + 1] - cp[k]), xp)
+
+        c = cm * 6
+
+        return xp.where(xp.abs(c) < eps, 0, c)
 
     @staticmethod
     @numba.jit(nopython=True, nogil=True, cache=True)
-    def _filter_periodic(s):  # pragma: no cover
+    def _filter_periodic_numba(s):  # pragma: no cover
         M = len(s)
         pole = -2 + np.sqrt(3)
 
         ndim = 1 if s.ndim == 1 else s.shape[1]
 
-        cp = np.zeros((M, ndim))
+        device = s.device if hasattr(s, "device") else None
+        cp = np.zeros((M, ndim), device=device)
         for k in range(M):
             cp[0] += s[(M - k) % M] * (pole**k)
         cp[0] *= 1 / (1 - (pole**M))
@@ -621,6 +973,38 @@ def _filter_periodic(s):  # pragma: no cover
         c = c.reshape(shape)
         return c
 
+    @staticmethod
+    def _filter_periodic_xp(s, xp):  # pragma: no cover
+        M = len(s)
+        pole = -2 + math.sqrt(3)
+
+        ndim = 1 if s.ndim == 1 else s.shape[1]
+
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(s.dtype, kind="integral") else s.dtype
+        device = s.device if hasattr(s, "device") else None
+        cp = xp.zeros((M, ndim), dtype=dtype, device=device)
+        for k in range(M):
+            cp = _item_add(cp, 0, s[(M - k) % M] * (pole**k), xp)
+        cp = _item_set(cp, 0, cp[0] * (1 / (1 - (pole**M))), xp)
+
+        for k in range(1, M):
+            cp = _item_set(cp, k, s[k] + pole * cp[k - 1], xp)
+
+        cm = xp.zeros((M, ndim), dtype=dtype, device=device)
+        for k in range(M):
+            cm = _item_add(cm, M - 1, (pole**k) * cp[k], xp)
+        cm = _item_set(cm, M - 1, cm[M - 1] * (pole / (1 - (pole**M))), xp)
+        cm = _item_add(cm, M - 1, cp[M - 1], xp)
+        cm = _item_set(cm, M - 1, cm[M - 1] * -pole, xp)
+
+        for k in range(M - 2, -1, -1):
+            cm = _item_set(cm, k, pole * (cm[k + 1] - cp[k]), xp)
+
+        c = cm * 6
+
+        eps = 1e-8
+        return xp.where(xp.abs(c) < eps, 0, c)
+
     def refinement_mask(self):
         order = int(self.support)
         mask = np.zeros(order + 1)
@@ -683,12 +1067,19 @@ def __repr__(self):
     def __eq__(self, other):
         return isinstance(other, type(self)) and other.M == self.M
 
-    def _func(self, t):
-        return self.__func(t, self.support / 2, self.M)
+    def _func_numba(self, t):
+        return self.__func_numba(t, self.support / 2, self.M)
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64, numba.float64, numba.float64)], nopython=True, cache=True)
-    def __func(t, half_support, M):  # pragma: no cover
+    @numba.vectorize(
+        [
+            numba.float32(numba.float32, numba.float32, numba.float32),
+            numba.float64(numba.float64, numba.float64, numba.float64),
+        ],
+        nopython=True,
+        cache=True,
+    )
+    def __func_numba(t, half_support, M):  # pragma: no cover
         t += half_support
 
         alpha = np.pi / M
@@ -704,12 +1095,19 @@ def __func(t, half_support, M):  # pragma: no cover
 
         return L * val
 
-    def _derivative_1(self, t):
-        return self.__derivative_1(t, self.support / 2, self.M)
+    def _derivative_1_numba(self, t):
+        return self.__derivative_1_numba(t, self.support / 2, self.M)
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64, numba.float64, numba.float64)], nopython=True, cache=True)
-    def __derivative_1(t, half_support, M):  # pragma: no cover
+    @numba.vectorize(
+        [
+            numba.float32(numba.float32, numba.float32, numba.float32),
+            numba.float64(numba.float64, numba.float64, numba.float64),
+        ],
+        nopython=True,
+        cache=True,
+    )
+    def __derivative_1_numba(t, half_support, M):  # pragma: no cover
         t += half_support
 
         alpha = np.pi / M
@@ -725,12 +1123,19 @@ def __derivative_1(t, half_support, M):  # pragma: no cover
 
         return L * val
 
-    def _derivative_2(self, t):
-        return self.__derivative_2(t, self.support / 2, self.M)
+    def _derivative_2_numba(self, t):
+        return self.__derivative_2_numba(t, self.support / 2, self.M)
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64, numba.float64, numba.float64)], nopython=True, cache=True)
-    def __derivative_2(t, half_support, M):  # pragma: no cover
+    @numba.vectorize(
+        [
+            numba.float32(numba.float32, numba.float32, numba.float32),
+            numba.float64(numba.float64, numba.float64, numba.float64),
+        ],
+        nopython=True,
+        cache=True,
+    )
+    def __derivative_2_numba(t, half_support, M):  # pragma: no cover
         t += half_support
 
         alpha = np.pi / M
@@ -750,20 +1155,98 @@ def __derivative_2(t, half_support, M):  # pragma: no cover
 
         return L * val
 
-    def filter_symmetric(self, s):
-        self.M = len(s)
+    def _func_xp(self, t, xp):  # pragma: no cover
+        half_support = self.support / 2
+        t = t + half_support
+
+        alpha = xp.pi / self.M
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(t.dtype, kind="integral") else t.dtype
+        device = t.device if hasattr(t, "device") else None
+        alpha = xp.asarray(alpha, dtype=dtype, device=device)
+        L = 1 / (4 * xp.sin(alpha) ** 2)
+
+        val = xp.where(
+            (t >= 0) & (t < 1),
+            2 * xp.sin(alpha * t) ** 2,
+            xp.where(
+                (t >= 1) & (t < 2),
+                xp.cos(2 * alpha * (t - 2)) + xp.cos(2 * alpha * (t - 1)) - 2 * xp.cos(2 * alpha),
+                xp.where(
+                    (t >= 2) & (t <= 3),
+                    2 * xp.sin(alpha * (t - 3)) ** 2,
+                    0.0,
+                ),
+            ),
+        )
+        return L * val
 
-        b0 = self._func(0)
-        b1 = self._func(1)
+    def _derivative_1_xp(self, t, xp):  # pragma: no cover
+        half_support = self.support / 2
+        t = t + half_support
+
+        alpha = xp.pi / self.M
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(t.dtype, kind="integral") else t.dtype
+        device = t.device if hasattr(t, "device") else None
+        alpha = xp.asarray(alpha, dtype=dtype, device=device)
+        L = 1 / (4 * xp.sin(alpha) ** 2)
+
+        val = xp.where(
+            (t >= 0) & (t <= 1),
+            4 * alpha * xp.sin(alpha * t) * xp.cos(alpha * t),
+            xp.where(
+                (t > 1) & (t <= 2),
+                2 * alpha * (xp.sin(2 * alpha * (2 - t)) + xp.sin(2 * alpha * (1 - t))),
+                xp.where(
+                    (t > 2) & (t <= 3),
+                    4 * alpha * xp.sin(alpha * (t - 3)) * xp.cos(alpha * (t - 3)),
+                    0.0,
+                ),
+            ),
+        )
+        return L * val
 
-        result = self._filter_symmetric(s.astype(float), self.M, b0, b1)
-        if s.ndim == 1:
-            return np.squeeze(result)
-        return result
+    def _derivative_2_xp(self, t, xp):  # pragma: no cover
+        half_support = self.support / 2
+        t = t + half_support
+
+        alpha = xp.pi / self.M
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(t.dtype, kind="integral") else t.dtype
+        device = t.device if hasattr(t, "device") else None
+        alpha = xp.asarray(alpha, dtype=dtype, device=device)
+        L = 1 / (4 * xp.sin(alpha) ** 2)
+
+        val = xp.where(
+            t < 0,
+            0.0,
+            xp.where(
+                (t > 0) & (t < 1),
+                4 * alpha**2 * (xp.cos(alpha * t) ** 2 - xp.sin(alpha * t) ** 2),
+                xp.where(
+                    (t > 1) & (t < 2),
+                    -4 * alpha**2 * (xp.cos(2 * alpha * (2 - t)) + xp.cos(2 * alpha * (1 - t))),
+                    xp.where(
+                        (t > 2) & (t < 3),
+                        4 * alpha**2 * (xp.cos(alpha * (t - 3)) ** 2 - xp.sin(alpha * (t - 3)) ** 2),
+                        xp.where(
+                            t > 3,
+                            0.0,
+                            xp.asarray(xp.nan, dtype=t.dtype),
+                        ),
+                    ),
+                ),
+            ),
+        )
+        return L * val
+
+    def _filter_symmetric_numba(self, s):
+        M = len(s)
+        b0 = self._func_numba(0)
+        b1 = self._func_numba(1)
+        return self.__filter_symmetric_numba(s, M, b0, b1)
 
     @staticmethod
     @numba.jit(nopython=True, nogil=True, cache=True)
-    def _filter_symmetric(s, M, b0, b1):  # pragma: no cover
+    def __filter_symmetric_numba(s, M, b0, b1):  # pragma: no cover
         ndim = 1 if s.ndim == 1 else s.shape[1]
         pole = (-b0 + np.sqrt(2 * b0 - 1)) / (1 - b0)
 
@@ -798,17 +1281,58 @@ def _filter_symmetric(s, M, b0, b1):  # pragma: no cover
         c = c.reshape(shape)
         return c
 
-    def filter_periodic(self, s):
-        self.M = len(s)
-        b0 = self._func(0)
-        result = self._filter_periodic(s.astype(float), self.M, b0)
-        if s.ndim == 1:
-            return np.squeeze(result)
-        return result
+    @staticmethod
+    def __b0_b1(M):
+        """
+        The values of the basis function at 0 and 1, computed with pure python
+        math so that the result can be compiled (e.g. with torch.compile).
+        """
+        alpha = math.pi / M
+        L = 1 / (4 * math.sin(alpha) ** 2)
+        b0 = L * (2 * math.cos(alpha) - 2 * math.cos(2 * alpha))
+        b1 = 2 * L * math.sin(alpha / 2) ** 2
+        return b0, b1
+
+    def _filter_symmetric_xp(self, s, xp):  # pragma: no cover
+        M = len(s)
+        b0, b1 = self.__b0_b1(self.M)
+        pole = (-b0 + math.sqrt(2 * b0 - 1)) / (1 - b0)
+
+        ndim = 1 if s.ndim == 1 else s.shape[1]
+
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(s.dtype, kind="integral") else s.dtype
+        device = s.device if hasattr(s, "device") else None
+        cp = xp.zeros((M, ndim), dtype=dtype, device=device)
+        eps = 1e-8
+        k0 = min(((2 * M) - 2, int(math.ceil(math.log(eps) / math.log(abs(pole))))))
+        for k in range(k0):
+            m = k % (2 * M - 2)
+            val = s[2 * M - 2 - m] if m >= M else s[m]
+            cp = _item_add(cp, 0, val * (pole**m), xp)
+        cp = _item_set(cp, 0, cp[0] * (1 / (1 - (pole ** (2 * M - 2)))), xp)
+
+        for k in range(1, M):
+            cp = _item_set(cp, k, s[k] + pole * cp[k - 1], xp)
+
+        cm = xp.zeros((M, ndim), dtype=dtype, device=device)
+        cm = _item_set(cm, M - 1, cp[M - 1] + (pole * cp[M - 2]), xp)
+        cm = _item_set(cm, M - 1, cm[M - 1] * (pole / ((pole * pole) - 1)), xp)
+
+        for k in range(M - 2, -1, -1):
+            cm = _item_set(cm, k, pole * (cm[k + 1] - cp[k]), xp)
+
+        c = cm / b1
+
+        return xp.where(xp.abs(c) < eps, 0, c)
+
+    def _filter_periodic_numba(self, s):
+        M = len(s)
+        b0 = self._func_numba(0)
+        return self.__filter_periodic_numba(s, M, b0)
 
     @staticmethod
     @numba.jit(nopython=True, nogil=True, cache=True)
-    def _filter_periodic(s, M, b0):  # pragma: no cover
+    def __filter_periodic_numba(s, M, b0):  # pragma: no cover
         ndim = 1 if s.ndim == 1 else s.shape[1]
 
         pole = (-b0 + np.sqrt(2 * b0 - 1)) / (1 - b0)
@@ -841,6 +1365,39 @@ def _filter_periodic(s, M, b0):  # pragma: no cover
         c = c.reshape(shape)
         return c
 
+    def _filter_periodic_xp(self, s, xp):  # pragma: no cover
+        M = len(s)
+        b0, _ = self.__b0_b1(self.M)
+        pole = (-b0 + math.sqrt(2 * b0 - 1)) / (1 - b0)
+
+        ndim = 1 if s.ndim == 1 else s.shape[1]
+
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(s.dtype, kind="integral") else s.dtype
+        device = s.device if hasattr(s, "device") else None
+        cp = xp.zeros((M, ndim), dtype=dtype, device=device)
+        cp = _item_set(cp, 0, s[0], xp)
+        for k in range(1, M):
+            cp = _item_add(cp, 0, s[k] * (pole ** (M - k)), xp)
+        cp = _item_set(cp, 0, cp[0] * (1 / (1 - (pole**M))), xp)
+
+        for k in range(1, M):
+            cp = _item_set(cp, k, s[k] + (pole * cp[k - 1]), xp)
+
+        cm = xp.zeros((M, ndim), dtype=dtype, device=device)
+        cm = _item_set(cm, M - 1, cp[M - 1], xp)
+        for k in range(M - 1):
+            cm = _item_add(cm, M - 1, cp[k] * (pole ** (k + 1)), xp)
+        cm = _item_set(cm, M - 1, cm[M - 1] * (1 / (1 - (pole**M))), xp)
+        cm = _item_set(cm, M - 1, cm[M - 1] * ((1 - pole) ** 2), xp)
+
+        for k in range(M - 2, -1, -1):
+            cm = _item_set(cm, k, (pole * cm[k + 1]) + (((1 - pole) ** 2) * cp[k]), xp)
+
+        c = cm
+
+        eps = 1e-8
+        return xp.where(xp.abs(c) < eps, 0, c)
+
     def refinement_mask(self):
         """
         TODO: This fails the tests. They are set to
@@ -914,8 +1471,8 @@ def __repr__(self):
         return "splinebox.basis_functions.CatmullRom()"
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _func(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _func_numba(t):  # pragma: no cover
         val = 0
         if np.abs(t) >= 0 and np.abs(t) <= 1:
             val = (3 / 2) * (np.abs(t) ** 3) - (5 / 2) * (np.abs(t) ** 2) + 1
@@ -924,8 +1481,8 @@ def _func(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _derivative_1(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _derivative_1_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t <= 1:
             val = t * (4.5 * t - 5)
@@ -938,8 +1495,8 @@ def _derivative_1(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def _derivative_2(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def _derivative_2_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t < 1:
             val = 9 * t - 5
@@ -954,12 +1511,81 @@ def _derivative_2(t):  # pragma: no cover
         return val
 
     @staticmethod
-    def filter_symmetric(s):
-        return s.astype(float)
+    def _func_xp(t, xp):  # pragma: no cover
+        abs_t = xp.abs(t)
+        return xp.where(
+            abs_t <= 1,
+            (3 / 2) * (abs_t**3) - (5 / 2) * (abs_t**2) + 1,
+            xp.where(
+                (abs_t > 1) & (abs_t <= 2),
+                (-1 / 2) * (abs_t**3) + (5 / 2) * (abs_t**2) - 4 * abs_t + 2,
+                0.0,
+            ),
+        )
 
     @staticmethod
-    def filter_periodic(s):
-        return s.astype(float)
+    def _derivative_1_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= 0) & (t <= 1),
+            t * (4.5 * t - 5),
+            xp.where(
+                (t >= -1) & (t < 0),
+                -t * (4.5 * t + 5),
+                xp.where(
+                    (t > 1) & (t <= 2),
+                    -1.5 * t * t + 5 * t - 4,
+                    xp.where(
+                        (t >= -2) & (t < -1),
+                        1.5 * t * t + 5 * t + 4,
+                        0.0,
+                    ),
+                ),
+            ),
+        )
+
+    @staticmethod
+    def _derivative_2_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            xp.abs(t) == 1,
+            xp.nan,
+            xp.where(
+                xp.abs(t) == 2,
+                xp.nan,
+                xp.where(
+                    (t >= 0) & (t < 1),
+                    9 * t - 5,
+                    xp.where(
+                        (t > -1) & (t < 0),
+                        -9 * t - 5,
+                        xp.where(
+                            (t > 1) & (t < 2),
+                            -3 * t + 5,
+                            xp.where(
+                                (t > -2) & (t < -1),
+                                3 * t + 5,
+                                0.0,
+                            ),
+                        ),
+                    ),
+                ),
+            ),
+        )
+
+    @staticmethod
+    def _filter_symmetric_numba(s):
+        return s
+
+    @staticmethod
+    def _filter_symmetric_xp(s, xp):
+        return s
+
+    @staticmethod
+    def _filter_periodic_numba(s):
+        return s
+
+    @staticmethod
+    def _filter_periodic_xp(s, xp):
+        return s
 
 
 class CubicHermite(BasisFunction):
@@ -1021,12 +1647,12 @@ def __str__(self):
     def __repr__(self):
         return "splinebox.basis_functions.CubicHermite()"
 
-    def _func(self, t):
-        return np.stack([self.h31(t), self.h32(t)], axis=-1)
+    def _func_numba(self, t):
+        return np.stack([self.h31_numba(t), self.h32_numba(t)], axis=-1)
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def h31(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def h31_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t <= 1:
             val = (1 + (2 * t)) * (t - 1) * (t - 1)
@@ -1035,8 +1661,8 @@ def h31(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def h32(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def h32_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t <= 1:
             val = t * (t - 1) * (t - 1)
@@ -1044,12 +1670,12 @@ def h32(t):  # pragma: no cover
             val = t * (t + 1) * (t + 1)
         return val
 
-    def _derivative_1(self, t):
-        return np.stack([self.h31prime(t), self.h32prime(t)], axis=-1)
+    def _derivative_1_numba(self, t):
+        return np.stack([self.h31prime_numba(t), self.h32prime_numba(t)], axis=-1)
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def h31prime(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def h31prime_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t <= 1:
             val = 6 * t * (t - 1)
@@ -1058,8 +1684,8 @@ def h31prime(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def h32prime(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def h32prime_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t <= 1:
             val = 3 * t * t - 4 * t + 1
@@ -1067,12 +1693,12 @@ def h32prime(t):  # pragma: no cover
             val = 3 * t * t + 4 * t + 1
         return val
 
-    def _derivative_2(self, t):
-        return np.stack([self.h31primeprime(t), self.h32primeprime(t)], axis=-1)
+    def _derivative_2_numba(self, t):
+        return np.stack([self.h31primeprime_numba(t), self.h32primeprime_numba(t)], axis=-1)
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def h31primeprime(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def h31primeprime_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t <= 1:
             val = 12 * t - 6
@@ -1081,8 +1707,8 @@ def h31primeprime(t):  # pragma: no cover
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64)], nopython=True, cache=True)
-    def h32primeprime(t):  # pragma: no cover
+    @numba.vectorize([numba.float32(numba.float32), numba.float64(numba.float64)], nopython=True, cache=True)
+    def h32primeprime_numba(t):  # pragma: no cover
         val = 0
         if t >= 0 and t <= 1:
             val = 6 * t - 4
@@ -1090,6 +1716,87 @@ def h32primeprime(t):  # pragma: no cover
             val = 6 * t + 4
         return val
 
+    def _func_xp(self, t, xp):  # pragma: no cover
+        return xp.stack([self.h31_xp(t, xp), self.h32_xp(t, xp)], axis=-1)
+
+    @staticmethod
+    def h31_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= 0) & (t <= 1),
+            (1 + (2 * t)) * (t - 1) * (t - 1),
+            xp.where(
+                (t < 0) & (t >= -1),
+                (1 - (2 * t)) * (t + 1) * (t + 1),
+                0.0,
+            ),
+        )
+
+    @staticmethod
+    def h32_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= 0) & (t <= 1),
+            t * (t - 1) * (t - 1),
+            xp.where(
+                (t < 0) & (t >= -1),
+                t * (t + 1) * (t + 1),
+                0.0,
+            ),
+        )
+
+    def _derivative_1_xp(self, t, xp):  # pragma: no cover
+        return xp.stack([self.h31prime_xp(t, xp), self.h32prime_xp(t, xp)], axis=-1)
+
+    @staticmethod
+    def h31prime_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= 0) & (t <= 1),
+            6 * t * (t - 1),
+            xp.where(
+                (t < 0) & (t >= -1),
+                -6 * t * (t + 1),
+                0.0,
+            ),
+        )
+
+    @staticmethod
+    def h32prime_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= 0) & (t <= 1),
+            3 * t * t - 4 * t + 1,
+            xp.where(
+                (t < 0) & (t >= -1),
+                3 * t * t + 4 * t + 1,
+                0.0,
+            ),
+        )
+
+    def _derivative_2_xp(self, t, xp):  # pragma: no cover
+        return xp.stack([self.h31primeprime_xp(t, xp), self.h32primeprime_xp(t, xp)], axis=-1)
+
+    @staticmethod
+    def h31primeprime_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= 0) & (t <= 1),
+            12 * t - 6,
+            xp.where(
+                (t < 0) & (t >= -1),
+                -12 * t - 6,
+                0.0,
+            ),
+        )
+
+    @staticmethod
+    def h32primeprime_xp(t, xp):  # pragma: no cover
+        return xp.where(
+            (t >= 0) & (t <= 1),
+            6 * t - 4,
+            xp.where(
+                (t < 0) & (t >= -1),
+                6 * t + 4,
+                0.0,
+            ),
+        )
+
     def h31_autocorrelation(self, i, j, M):  # pragma: no cover
         """
         Derived by V. Uhlmann during her PhD.
@@ -1202,12 +1909,20 @@ def h3_periodic_crosscorrelation(self, n, M):  # pragma: no cover
         return val
 
     @staticmethod
-    def filter_symmetric(s):
-        return s.astype(float)
+    def _filter_symmetric_numba(s):
+        return s
 
     @staticmethod
-    def filter_periodic(s):
-        return s.astype(float)
+    def _filter_symmetric_xp(s, xp):
+        return s
+
+    @staticmethod
+    def _filter_periodic_numba(s):
+        return s
+
+    @staticmethod
+    def _filter_periodic_xp(s, xp):
+        return s
 
 
 class ExponentialHermite(BasisFunction):
@@ -1273,12 +1988,16 @@ def __repr__(self):
     def __eq__(self, other):
         return isinstance(other, type(self)) and other.M == self.M
 
-    def _func(self, t):
-        return np.stack([self._he31(t, self.M), self._he32(t, self.M)], axis=-1)
+    def _func_numba(self, t):
+        return np.stack([self._he31_numba(t, self.M), self._he32_numba(t, self.M)], axis=-1)
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64, numba.float64)], nopython=True, cache=True)
-    def _he31(t, M):  # pragma: no cover
+    @numba.vectorize(
+        [numba.float32(numba.float32, numba.float32), numba.float64(numba.float64, numba.float64)],
+        nopython=True,
+        cache=True,
+    )
+    def _he31_numba(t, M):  # pragma: no cover
         def _g1(t, M):
             val = 0
             if t >= 0 and t <= 1:
@@ -1296,8 +2015,12 @@ def _g1(t, M):
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64, numba.float64)], nopython=True, cache=True)
-    def _he32(t, M):  # pragma: no cover
+    @numba.vectorize(
+        [numba.float32(numba.float32, numba.float32), numba.float64(numba.float64, numba.float64)],
+        nopython=True,
+        cache=True,
+    )
+    def _he32_numba(t, M):  # pragma: no cover
         def _g2(t, M):
             val = 0
             if t >= 0 and t <= 1:
@@ -1315,12 +2038,16 @@ def _g2(t, M):
         val = _g2(t, M) if t >= 0 else -1 * _g2(-t, M)
         return val
 
-    def _derivative_1(self, t):
-        return np.stack([self._he31prime(t, self.M), self._he32prime(t, self.M)], axis=-1)
+    def _derivative_1_numba(self, t):
+        return np.stack([self._he31prime_numba(t, self.M), self._he32prime_numba(t, self.M)], axis=-1)
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64, numba.float64)], nopython=True, cache=True)
-    def _he31prime(t, M):  # pragma: no cover
+    @numba.vectorize(
+        [numba.float32(numba.float32, numba.float32), numba.float64(numba.float64, numba.float64)],
+        nopython=True,
+        cache=True,
+    )
+    def _he31prime_numba(t, M):  # pragma: no cover
         def _g1prime(t, M):
             val = 0
             if t >= 0 and t <= 1:
@@ -1334,8 +2061,12 @@ def _g1prime(t, M):
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64, numba.float64)], nopython=True, cache=True)
-    def _he32prime(t, M):  # pragma: no cover
+    @numba.vectorize(
+        [numba.float32(numba.float32, numba.float32), numba.float64(numba.float64, numba.float64)],
+        nopython=True,
+        cache=True,
+    )
+    def _he32prime_numba(t, M):  # pragma: no cover
         def _g2prime(t, M):
             val = 0
             if t >= 0 and t <= 1:
@@ -1352,12 +2083,16 @@ def _g2prime(t, M):
         val = _g2prime(t, M) if t >= 0 else _g2prime(-t, M)
         return val
 
-    def _derivative_2(self, t):
-        return np.stack([self._he31primeprime(t, self.M), self._he32primeprime(t, self.M)], axis=-1)
+    def _derivative_2_numba(self, t):
+        return np.stack([self._he31primeprime_numba(t, self.M), self._he32primeprime_numba(t, self.M)], axis=-1)
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64, numba.float64)], nopython=True, cache=True)
-    def _he31primeprime(t, M):  # pragma: no cover
+    @numba.vectorize(
+        [numba.float32(numba.float32, numba.float32), numba.float64(numba.float64, numba.float64)],
+        nopython=True,
+        cache=True,
+    )
+    def _he31primeprime_numba(t, M):  # pragma: no cover
         def _g1primeprime(t, M):
             val = 0
             if t >= 0 and t <= 1:
@@ -1371,8 +2106,12 @@ def _g1primeprime(t, M):
         return val
 
     @staticmethod
-    @numba.vectorize([numba.float64(numba.float64, numba.float64)], nopython=True, cache=True)
-    def _he32primeprime(t, M):  # pragma: no cover
+    @numba.vectorize(
+        [numba.float32(numba.float32, numba.float32), numba.float64(numba.float64, numba.float64)],
+        nopython=True,
+        cache=True,
+    )
+    def _he32primeprime_numba(t, M):  # pragma: no cover
         def _g2primeprime(t, M):
             val = 0
             if t >= 0 and t <= 1:
@@ -1387,13 +2126,114 @@ def _g2primeprime(t, M):
         val = _g2primeprime(t, M) if t >= 0 else -1 * _g2primeprime(-t, M)
         return val
 
+    def _func_xp(self, t, xp):  # pragma: no cover
+        return xp.stack([self._he31_xp(t, self.M, xp), self._he32_xp(t, self.M, xp)], axis=-1)
+
+    @staticmethod
+    def _he31_xp(t, M, xp):  # pragma: no cover
+        abs_t = xp.abs(t)
+        alpha = xp.pi / M
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(t.dtype, kind="integral") else t.dtype
+        device = t.device if hasattr(t, "device") else None
+        alpha = xp.asarray(alpha, dtype=dtype, device=device)
+        denom = (alpha * xp.cos(alpha)) - xp.sin(alpha)
+        num = (
+            (0.5 * (2 * alpha * xp.cos(alpha) - xp.sin(alpha)))
+            - (alpha * xp.cos(alpha) * abs_t)
+            - (0.5 * xp.sin(alpha - (2 * alpha * abs_t)))
+        )
+        return xp.where((abs_t >= 0) & (abs_t <= 1), num / denom, 0.0)
+
+    @staticmethod
+    def _he32_xp(t, M, xp):  # pragma: no cover
+        abs_t = xp.abs(t)
+        alpha = xp.pi / M
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(t.dtype, kind="integral") else t.dtype
+        device = t.device if hasattr(t, "device") else None
+        alpha = xp.asarray(alpha, dtype=dtype, device=device)
+        denom = ((alpha * xp.cos(alpha)) - xp.sin(alpha)) * 8 * alpha * xp.sin(alpha)
+        num = (
+            -((2 * alpha * xp.cos(2 * alpha)) - xp.sin(2 * alpha))
+            - (4 * alpha * xp.sin(alpha) * xp.sin(alpha) * abs_t)
+            - (2 * xp.sin(alpha) * xp.cos(2 * alpha * (abs_t - 0.5)))
+            + (2 * alpha * xp.cos(2 * alpha * (abs_t - 1)))
+        )
+        val = xp.where((abs_t >= 0) & (abs_t <= 1), num / denom, 0.0)
+        return xp.where(t >= 0, val, -val)
+
+    def _derivative_1_xp(self, t, xp):  # pragma: no cover
+        return xp.stack([self._he31prime_xp(t, self.M, xp), self._he32prime_xp(t, self.M, xp)], axis=-1)
+
+    @staticmethod
+    def _he31prime_xp(t, M, xp):  # pragma: no cover
+        abs_t = xp.abs(t)
+        alpha = xp.pi / M
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(t.dtype, kind="integral") else t.dtype
+        device = t.device if hasattr(t, "device") else None
+        alpha = xp.asarray(alpha, dtype=dtype, device=device)
+        denom = (alpha * xp.cos(alpha)) - xp.sin(alpha)
+        num = -(alpha * xp.cos(alpha)) + (alpha * xp.cos(alpha - (2 * alpha * abs_t)))
+        val = xp.where((abs_t >= 0) & (abs_t <= 1), num / denom, 0.0)
+        return xp.where(t >= 0, val, -val)
+
+    @staticmethod
+    def _he32prime_xp(t, M, xp):  # pragma: no cover
+        abs_t = xp.abs(t)
+        alpha = xp.pi / M
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(t.dtype, kind="integral") else t.dtype
+        device = t.device if hasattr(t, "device") else None
+        alpha = xp.asarray(alpha, dtype=dtype, device=device)
+        denom = ((alpha * xp.cos(alpha)) - xp.sin(alpha)) * 8 * alpha * xp.sin(alpha)
+        num = (
+            -(4 * alpha * xp.sin(alpha) * xp.sin(alpha))
+            + (4 * alpha * xp.sin(alpha) * xp.sin(2 * alpha * (abs_t - 0.5)))
+            - (4 * alpha**2 * xp.sin(2 * alpha * (abs_t - 1)))
+        )
+        return xp.where((abs_t >= 0) & (abs_t <= 1), num / denom, 0.0)
+
+    def _derivative_2_xp(self, t, xp):  # pragma: no cover
+        return xp.stack([self._he31primeprime_xp(t, self.M, xp), self._he32primeprime_xp(t, self.M, xp)], axis=-1)
+
+    @staticmethod
+    def _he31primeprime_xp(t, M, xp):  # pragma: no cover
+        abs_t = xp.abs(t)
+        alpha = xp.pi / M
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(t.dtype, kind="integral") else t.dtype
+        device = t.device if hasattr(t, "device") else None
+        alpha = xp.asarray(alpha, dtype=dtype, device=device)
+        denom = (alpha * xp.cos(alpha)) - xp.sin(alpha)
+        num = 2 * alpha**2 * xp.sin(alpha - (2 * alpha * abs_t))
+        return xp.where((abs_t >= 0) & (abs_t <= 1), num / denom, 0.0)
+
+    @staticmethod
+    def _he32primeprime_xp(t, M, xp):  # pragma: no cover
+        abs_t = xp.abs(t)
+        alpha = xp.pi / M
+        dtype = xp.asarray(1.0).dtype if xp.isdtype(t.dtype, kind="integral") else t.dtype
+        device = t.device if hasattr(t, "device") else None
+        alpha = xp.asarray(alpha, dtype=dtype, device=device)
+        denom = ((alpha * xp.cos(alpha)) - xp.sin(alpha)) * 8 * alpha * xp.sin(alpha)
+        num = +(8 * alpha**2 * xp.sin(alpha) * xp.cos(2 * alpha * (abs_t - 0.5))) - (
+            8 * alpha**3 * xp.cos(2 * alpha * (abs_t - 1))
+        )
+        val = xp.where((abs_t >= 0) & (abs_t <= 1), num / denom, 0.0)
+        return xp.where(t >= 0, val, -val)
+
+    @staticmethod
+    def _filter_symmetric_numba(s):
+        return s
+
+    @staticmethod
+    def _filter_symmetric_xp(s, xp):
+        return s
+
     @staticmethod
-    def filter_symmetric(s):
-        return s.astype(float)
+    def _filter_periodic_numba(s):
+        return s
 
     @staticmethod
-    def filter_periodic(s):
-        return s.astype(float)
+    def _filter_periodic_xp(s, xp):
+        return s
 
 
 def _multinomial(
diff --git a/src/splinebox/spline_curves.py b/src/splinebox/spline_curves.py
index 9c770dc..02c5be0 100644
--- a/src/splinebox/spline_curves.py
+++ b/src/splinebox/spline_curves.py
@@ -7,13 +7,17 @@
 import itertools
 import json
 import math
+import numbers
+import re
 import warnings
 
+import array_api_compat
 import numba
 import numpy as np
 import scipy.integrate
 
 import splinebox.basis_functions
+from splinebox.basis_functions import _item_set
 
 GAUSS_LEGENDRE_QUADRATURE_POINTS = np.array(
     [
@@ -80,8 +84,12 @@ def padding_function(knots, pad_length):
            [3, 3],
            [3, 3]])
     """
+    xp = array_api_compat.array_namespace(knots)
+
     # Add constant padding to the ends
-    return np.pad(knots, ((pad_length, pad_length), (0, 0)), mode="edge")
+    before = xp.tile(knots[0], (pad_length, 1))
+    after = xp.tile(knots[-1], (pad_length, 1))
+    return xp.concat((before, knots, after), axis=0)
 
 
 class Spline:
@@ -154,7 +162,6 @@ class Spline:
     _wrong_array_size_msg = (
         "It looks like control_points is neither a 1 nor a 2D array. I don't know how to handle this yet."
     )
-    _no_control_points_msg = "This spline object doesn't have any control points yet."
     _unimplemented_msg = "This function is not implemented."
 
     def __init__(
@@ -188,14 +195,16 @@ def __init__(
         self._cached_segments = None
         self._cached_segment_lengths = None
 
-    def _check_control_points(self):
+    def _get_control_points_xp_and_device(self):
         """
         Most methods require control points to be set before they
         can be used. This helper function checks if the control points have been
-        set.
+        set and return the their array namespace.
         """
         if self.control_points is None:
-            raise RuntimeError(self._no_control_points_msg)
+            raise RuntimeError("This spline object doesn't have any control points yet.")
+        else:
+            return array_api_compat.array_namespace(self.control_points), array_api_compat.device(self.control_points)
 
     def __str__(self):
         closed_str = "closed" if self.closed else "open"
@@ -205,14 +214,37 @@ def __str__(self):
             return f"{closed_str} {self.ndim}D {self.basis_function} spline with {self.M} knots"
 
     def __repr__(self):
-        return f"splinebox.spline_curves.Spline(M={repr(self.M)}, basis_function={repr(self.basis_function)}, closed={repr(self.closed)}, control_points=np.{repr(self.control_points)})"
+        xp, device = self._get_control_points_xp_and_device()
+        control_points_repr = repr(self.control_points)
+        control_points_repr = (
+            re.split(r",\s*dtype\b", control_points_repr)[0]
+            if "dtype" in control_points_repr
+            else control_points_repr[:-1]
+        )
+        dtype_repr = f"{xp.__name__}.{repr(self.control_points.dtype)}"
+
+        if array_api_compat.is_numpy_namespace(xp):
+            control_points_repr = f"np.{control_points_repr}, dtype={dtype_repr})"
+        elif array_api_compat.is_jax_namespace(xp):
+            device_repr = f"jax.devices('{device.platform}')[{device.id}]"
+            control_points_repr = control_points_repr.replace("Array", "array")
+            control_points_repr = f"{xp.__name__}.{control_points_repr}, dtype={dtype_repr}, device={device_repr})"
+        else:
+            device_repr = f"{xp.__name__}.{repr(device)}"
+            control_points_repr = f"{xp.__name__}.{control_points_repr}, dtype={dtype_repr}, device={device_repr})"
+
+        return f"splinebox.spline_curves.Spline(M={repr(self.M)}, basis_function={repr(self.basis_function)}, closed={repr(self.closed)}, control_points={control_points_repr})"
 
     def __eq__(self, other):
+        xp, device = self._get_control_points_xp_and_device()
         return (
             self.M == other.M
             and self.basis_function == other.basis_function
             and self.closed == other.closed
-            and np.all(self.control_points == other.control_points)
+            and self.closed == other.closed
+            and xp.all(self.control_points == other.control_points)
+            if self.control_points is not None and other.control_points is not None
+            else self.control_points is other.control_points
         )
 
     @property
@@ -231,7 +263,8 @@ def control_points(self):
     @control_points.setter
     def control_points(self, values):
         if values is not None:
-            values = np.array(values)
+            if not array_api_compat.is_array_api_obj(values):
+                values = np.array(values)
             if values.ndim != 2:
                 raise ValueError(
                     "The control point array has to be 2D. The first dimension encodes the control points and the second corresponds to the dimensionality of the codomain of the spline."
@@ -261,15 +294,20 @@ def knots(self):
         ValueError
             If the number of knots doesn't match :attr:`~splinebox.spline_curves.Spline.M`.
         """
+        xp, device = self._get_control_points_xp_and_device()
         if self.padding_function is None and not self.closed:
-            t = np.arange(-self.pad, self.M + self.pad)
+            t = xp.arange(-self.pad, self.M + self.pad, device=device)
         else:
-            t = np.arange(self.M)
+            t = xp.arange(self.M, device=device)
         return self(t)
 
     @knots.setter
     def knots(self, values):
-        knots = np.array(values)
+        if not array_api_compat.is_array_api_obj(values):
+            raise ValueError(
+                "The knots have to be provided as an array API compatible array, e.g. numpy array, torch tensor, ..."
+            )
+        knots = values
         if knots.ndim != 2:
             raise ValueError(
                 "The knots array has to be 2D. The first dimension encodes the knots and the second corresponds to the dimensionality of the codomain of the spline."
@@ -415,16 +453,23 @@ def integration_segment_size(self, value):
     @property
     def _segments(self):
         if self._cached_segments is None:
+            xp, device = self._get_control_points_xp_and_device()
             if self.closed:
-                self._cached_segments = np.arange(
-                    0, self.M + self.integration_segment_size / 10, self.integration_segment_size
+                self._cached_segments = xp.arange(
+                    0,
+                    self.M + self.integration_segment_size / 10,
+                    self.integration_segment_size,
+                    device=device,
+                    dtype=xp.result_type(self.control_points, xp.float32),
                 )
             else:
                 self._cached_segments = (
-                    np.arange(
+                    xp.arange(
                         0,
                         self.M - 1 + 2 * (self.pad + self.half_support) + self.integration_segment_size / 10,
                         self.integration_segment_size,
+                        device=device,
+                        dtype=xp.result_type(self.control_points, xp.float32),
                     )
                     - self.pad
                     - self.half_support
@@ -434,14 +479,17 @@ def _segments(self):
     @property
     def _segment_lengths(self):
         if self._cached_segment_lengths is None:
-            tvals = np.add.outer(
-                (self._segments[:-1] + self._segments[1:]) / 2,
-                GAUSS_LEGENDRE_QUADRATURE_POINTS / (2 / self.integration_segment_size),
+            xp, device = self._get_control_points_xp_and_device()
+            tvals = xp.expand_dims((self._segments[:-1] + self._segments[1:]) / 2, axis=1) + xp.expand_dims(
+                xp.asarray(GAUSS_LEGENDRE_QUADRATURE_POINTS, device=device) / (2 / self.integration_segment_size),
+                axis=0,
             )
             function_values = self(tvals.flatten(), derivative=1).reshape(*tvals.shape, -1)
-            function_values = np.linalg.norm(function_values, axis=-1)
+            function_values = xp.linalg.vector_norm(function_values, axis=-1)
             self._cached_segment_lengths = (
-                function_values @ GAUSS_LEGENDRE_QUADRATURE_WEIGHTS / (2 / self.integration_segment_size)
+                function_values
+                @ xp.asarray(GAUSS_LEGENDRE_QUADRATURE_WEIGHTS, device=device)
+                / (2 / self.integration_segment_size)
             )
         return self._cached_segment_lengths
 
@@ -796,7 +844,7 @@ def draw(self, x, y):
         >>> plt.plot(vals[:, 0], vals[:, 1])  # doctest: +SKIP
         >>> plt.show()  # doctest: +SKIP
         """
-        self._check_control_points()
+        xp, device = self._get_control_points_xp_and_device()
 
         if self.ndim != 2:
             raise RuntimeError("draw() can only be used with 2D curves")
@@ -804,7 +852,7 @@ def draw(self, x, y):
         if not self.closed:
             raise RuntimeError("draw() can only be used with closed curves.")
 
-        xx, yy = np.meshgrid(x, y)
+        xx, yy = xp.meshgrid(x, y)
         result = self.is_inside(xx.flatten(), yy.flatten())
         return result.reshape(xx.shape)
 
@@ -822,22 +870,22 @@ def dtheta(self, t):
         .. math::
             \frac{d \theta}{dt} = \frac{1}{r^2} \left( x\frac{dy}{dt} - y\frac{dx}{dt} \right) \text{, where } r^2 = x^2 + y^2
         """
-        self._check_control_points()
+        xp, device = self._get_control_points_xp_and_device()
         if self.ndim != 2:
             raise RuntimeError("dtheta() is only defined for 2D curves.")
-        t, single_value = self._convert_to_array(t)
+        t, single_value = self._convert_to_array(t, xp, device)
         r = self(t)
         dr = self(t, derivative=1)
         if r.ndim == 1:
-            r = r[np.newaxis, :]
+            r = r[xp.newaxis, :]
             dr = dr[np.newaxis, :]
-        r2 = np.linalg.norm(r, axis=1) ** 2
+        r2 = xp.linalg.norm(r, axis=1) ** 2
         val = (1.0 / r2) * (r[:, 0] * dr[:, 1] - r[:, 1] * dr[:, 0])
         # Happens if the spline is not differentiable in this location
         # or when the point is at the origin
         # The solution to return 0 is a bit of a hack but works in practice with
         # the numerical integration
-        val = np.nan_to_num(val)
+        val = xp.nan_to_num(val)
         if single_value:
             val = val[0]
         return val
@@ -891,29 +939,46 @@ def is_inside(self, x, y):
         """
         if not self.closed:
             raise RuntimeError("isInside() can only be used with closed curves.")
-        self._check_control_points()
+
+        xp, device = self._get_control_points_xp_and_device()
+
         if self.ndim != 2:
             raise RuntimeError("isInside() can only be used with 2D curves.")
 
-        if isinstance(x, (float, int)):
-            x = np.array([x])
-        if isinstance(y, (float, int)):
-            y = np.array([y])
+        if isinstance(x, numbers.Number):
+            x = xp.asarray([x], device=device)
+        if isinstance(y, numbers.Number):
+            y = xp.asarray([y], device=device)
+
+        if array_api_compat.array_namespace(x) is not xp or array_api_compat.array_namespace(y) is not xp:
+            raise TypeError(
+                f"The control points are a {xp.__name__.split('.')[-1]} array but "
+                f"x is a {array_api_compat.array_namespace(x).__name__.split('.')[-1]} array "
+                f"and y is a {array_api_compat.array_namespace(y).__name__.split('.')[-1]} array."
+            )
+
         if not np.allclose(x.shape, y.shape):
             raise ValueError("x and y need to have the same shape")
 
-        results = np.zeros(x.shape)
-        for coord, point in enumerate(np.stack([x, y], axis=-1)):
+        results = xp.zeros(x.shape, device=device)
+        for coord, point in enumerate(xp.stack([x, y], axis=-1)):
             spline_copy = self.copy()
             spline_copy.translate(-point)
             winding_number = scipy.integrate.quad(spline_copy.dtheta, 0, spline_copy.M)[0]
             if np.abs(np.abs(winding_number) - 2 * np.pi) < 1e-6:
-                results[coord] = 1
+                results = _item_set(results, coord, 1, xp)
             elif np.abs(winding_number) > 1e-6:
-                results[coord] = 0.5
+                results = _item_set(results, coord, 0.5, xp)
             else:
-                results[coord] = 0
-        results = np.squeeze(results)
+                results = _item_set(results, coord, 0, xp)
+
+        # Squeeze out all dimension that are 1
+        axes = tuple(i for i, size in enumerate(results.shape) if size == 1)
+        if axes:
+            results = xp.squeeze(
+                results,
+            )
+
         if results.ndim == 0:
             return float(results)
         return results
@@ -956,7 +1021,13 @@ def fit(self, points, boundary_condition="free"):
                [ 8.   ,  2.811],
                [10.333, -4.283]])
         """
-        if len(points) < self.M + 2 * self.pad:
+        if not array_api_compat.is_array_api_obj(points):
+            raise TypeError("points has to be an array api compatible object, e.g. numpy array, torch tensor, ...")
+        xp = array_api_compat.array_namespace(points)
+        device = array_api_compat.device(points)
+
+        n_points = len(points)
+        if n_points < self.M + 2 * self.pad:
             raise RuntimeError(
                 f"You provided too few points. For the fit to have a unique solution you need to provide at least as many points as your spline has control points (including padding). This spline has {self.M}+{2 * self.pad}(padding) control points but you provided {len(points)} points. Consider providing more points or reducing the number of knots M. If the number of points is equal to M, consider using `spline.knots = points` instead of fitting."
             )
@@ -965,7 +1036,8 @@ def fit(self, points, boundary_condition="free"):
                 "The points should be provided as a 2D array, with the first dimension encoding the different data points and the second dimension encoding the dimensionality of the data points."
             )
 
-        n_points = len(points)
+        # We do the fit in numpy to be able to exploit scipy's sparse arrays
+        points = np.array(points)
         n_control_points = self.M if self.closed else self.M + 2 * self.pad
 
         t = np.linspace(0, self.M, n_points + 1)[:-1] if self.closed else np.linspace(0, self.M - 1, n_points)
@@ -974,7 +1046,7 @@ def fit(self, points, boundary_condition="free"):
         shift = np.arange(-bound + 1, bound + 1)
         tval = np.empty((len(t), len(shift)), dtype=float)
         indices = np.empty((len(t), len(shift)), dtype=int)
-        self._compute_tval_and_indices(t, shift, self.closed, self.M, self.pad, tval, indices)
+        self._compute_tval_and_indices_numba(t, shift, self.closed, self.M, self.pad, tval, indices)
 
         row = np.repeat(np.arange(indices.shape[0]), indices.shape[1])
         col = indices.flatten()
@@ -989,8 +1061,7 @@ def fit(self, points, boundary_condition="free"):
                 data = data[mask]
             basis_function_values = scipy.sparse.csr_array((data, (row, col)), shape=(n_points, n_control_points))
 
-            if self.control_points is None:
-                self.control_points = np.empty((n_control_points, points.shape[1]))
+            self.control_points = np.empty((n_control_points, points.shape[1]))
             for i in range(self.ndim):
                 self.control_points[:, i] = scipy.sparse.linalg.lsqr(basis_function_values, points[:, i])[0]
         elif boundary_condition in ("clamped", "natural"):
@@ -1018,10 +1089,9 @@ def fit(self, points, boundary_condition="free"):
             )
             basis_function_values = scipy.sparse.csr_array((data, (row, col)), shape=(n_points, n_control_points - 2))
 
-            if self.control_points is None:
-                if points.ndim == 1:
-                    points = points[:, np.newaxis]
-                self.control_points = np.empty((n_control_points, points.shape[1]))
+            if points.ndim == 1:
+                points = points[:, np.newaxis]
+            self.control_points = np.empty((n_control_points, points.shape[1]))
             for i in range(self.ndim):
 
                 self.control_points[1:-1, i] = scipy.sparse.linalg.lsqr(basis_function_values, points[:, i])[0]
@@ -1038,13 +1108,21 @@ def fit(self, points, boundary_condition="free"):
         else:
             raise ValueError(f"Unknown boundary_conditions {boundary_condition}")
 
+        # Convert to the correct array type
+        if not array_api_compat.is_numpy_namespace(xp):
+            self.control_points = xp.asarray(self.control_points, device=device)
+
     def _gauss_legendre_quadrature(self, a, b):
-        tvals = np.outer((b - a) / 2, GAUSS_LEGENDRE_QUADRATURE_POINTS) + (a + b)[:, np.newaxis] / 2
+        xp, device = self._get_control_points_xp_and_device()
+        tvals = (
+            xp.outer((b - a) / 2, xp.asarray(GAUSS_LEGENDRE_QUADRATURE_POINTS, device=device))
+            + (a + b)[:, np.newaxis] / 2
+        )
         function_values = self(tvals.flatten(), derivative=1)
         function_values = function_values.reshape(*tvals.shape, -1)
-        function_values = np.linalg.norm(function_values, axis=-1)
-        result = function_values @ GAUSS_LEGENDRE_QUADRATURE_WEIGHTS * (b - a) / 2
-        result[a == b] = 0
+        function_values = xp.linalg.vector_norm(function_values, axis=-1)
+        result = function_values @ xp.asarray(GAUSS_LEGENDRE_QUADRATURE_WEIGHTS, device=device) * (b - a) / 2
+        result = xp.where(a == b, 0, result)
         return result
 
     def arc_length(self, stop=None, start=0):
@@ -1090,54 +1168,59 @@ def arc_length(self, stop=None, start=0):
         >>> spline.arc_length(np.arange(M - 1), np.arange(1, M))
         array([4.982, 5.288, 3.09 , 3.675])
         """
-        self._check_control_points()
+        xp, device = self._get_control_points_xp_and_device()
 
         if stop is None:
             stop = self.M if self.closed else self.M - 1
 
-        start, start_single_value = self._convert_to_array(start)
-        stop, stop_single_value = self._convert_to_array(stop)
+        start, start_single_value = self._convert_to_array(start, xp, device)
+        stop, stop_single_value = self._convert_to_array(stop, xp, device)
 
         single_value = start_single_value and stop_single_value
 
         if not single_value:
             if start_single_value:
-                start = np.repeat(start, len(stop))
+                start = xp.repeat(start, len(stop))
             elif stop_single_value:
-                stop = np.repeat(stop, len(start))
+                stop = xp.repeat(stop, len(start))
             else:
                 if len(stop) != len(start):
                     raise ValueError(
                         "If you provide array like objects for start and stop, they need to have the same length."
                     )
-        bounds = np.stack([start, stop], axis=-1)
-        bounds = np.sort(bounds, axis=-1)
+        bounds = xp.stack([start, stop], axis=-1)
+        bounds = xp.sort(bounds, axis=-1)
 
-        start_segment_indices = np.searchsorted(self._segments, bounds[:, 0], side="left")
-        stop_segment_indices = np.searchsorted(self._segments, bounds[:, 1], side="right") - 1
+        start_segment_indices = xp.searchsorted(self._segments, bounds[:, 0], side="left")
+        stop_segment_indices = xp.searchsorted(self._segments, bounds[:, 1], side="right") - 1
         initial_segment_lengths = self._gauss_legendre_quadrature(bounds[:, 0], self._segments[start_segment_indices])
         final_segment_lengths = self._gauss_legendre_quadrature(self._segments[stop_segment_indices], bounds[:, 1])
 
-        arc_lengths = np.empty(len(bounds))
+        arc_lengths = xp.empty(len(bounds), dtype=xp.result_type(self.control_points, bounds))
         for i, (initial_segment_length, final_segment_length) in enumerate(
             zip(initial_segment_lengths, final_segment_lengths)
         ):
             if start_segment_indices[i] > stop_segment_indices[i]:
                 # Happens when both bounds are within the same segment
-                arc_lengths[i] = self._gauss_legendre_quadrature(
-                    np.array([bounds[i, 0]]),
-                    np.array([bounds[i, 1]]),
-                )[0]
+                arc_lengths = _item_set(
+                    arc_lengths,
+                    i,
+                    self._gauss_legendre_quadrature(
+                        xp.asarray([bounds[i, 0]], device=device),
+                        xp.asarray([bounds[i, 1]], device=device),
+                    )[0],
+                    xp,
+                )
             else:
                 segment_sum = math.fsum(self._segment_lengths[start_segment_indices[i] : stop_segment_indices[i]])
-                arc_lengths[i] = segment_sum + initial_segment_length + final_segment_length
+                arc_lengths = _item_set(arc_lengths, i, segment_sum + initial_segment_length + final_segment_length, xp)
 
         if single_value:
             arc_lengths = arc_lengths[0]
         return arc_lengths
 
     def _length_to_parameter_recursion(
-        self, s, current_value, lower_bound, upper_bound, intermediate_results=None, atol=1e-4
+        self, s, current_value, lower_bound, upper_bound, xp, intermediate_results=None, atol=1e-4
     ):
         """
         Convert the given arc length s on the curve to a value in parameter space.
@@ -1169,16 +1252,15 @@ def _length_to_parameter_recursion(
         midpoint_length = current_value + self.arc_length(lower_bound, midpoint)
         if intermediate_results is not None:
             intermediate_results.append((midpoint, midpoint_length))
-
-        if np.isclose(s, midpoint_length, atol=atol, rtol=0):
+        if xp.isclose(s, midpoint_length, atol=atol, rtol=0):
             return midpoint
         elif s < midpoint_length:
             return self._length_to_parameter_recursion(
-                s, current_value, lower_bound, midpoint, intermediate_results, atol
+                s, current_value, lower_bound, midpoint, xp, intermediate_results, atol
             )
         else:
             return self._length_to_parameter_recursion(
-                s, midpoint_length, midpoint, upper_bound, intermediate_results, atol
+                s, midpoint_length, midpoint, upper_bound, xp, intermediate_results, atol
             )
 
     def arc_length_to_parameter(self, s, atol=1e-4):
@@ -1209,14 +1291,15 @@ def arc_length_to_parameter(self, s, atol=1e-4):
         >>> spline.arc_length(0, 2.12)  # doctest: +NUMBER
         2.2
         """
-        self._check_control_points()
-        s, single_value = self._convert_to_array(s)
-        sort_indices = np.argsort(s)
-        results = np.zeros_like(s, dtype=float)
+        xp, device = self._get_control_points_xp_and_device()
+        s, single_value = self._convert_to_array(s, xp, device)
+        sort_indices = xp.argsort(s)
+        results = xp.empty(s.shape, dtype=xp.result_type(self.control_points, s))
+        s = xp.astype(s, results.dtype)
 
-        current_value = 0
-        lower_bound = 0
-        upper_bound = self.M if self.closed else self.M - 1
+        current_value = xp.asarray(0, dtype=s.dtype)
+        lower_bound = xp.asarray(0, dtype=s.dtype)
+        upper_bound = xp.asarray(self.M if self.closed else self.M - 1, dtype=s.dtype)
         intermediate_results = []
 
         total_length = self.arc_length()
@@ -1224,7 +1307,7 @@ def arc_length_to_parameter(self, s, atol=1e-4):
         def upper_bound_key_func(x):
             diff = x[1] - s[i]
             if diff < 0:
-                return np.inf
+                return xp.asarray(xp.inf, dtype=s.dtype)
             return diff
 
         for i in sort_indices:
@@ -1235,10 +1318,15 @@ def upper_bound_key_func(x):
             if len(intermediate_results) > 0:
                 upper_bound, upper_bound_len = min(intermediate_results, key=upper_bound_key_func)
                 if s[i] > upper_bound_len:
-                    upper_bound = self.M if self.closed else self.M - 1
-
-            results[i] = self._length_to_parameter_recursion(
-                s[i], current_value, lower_bound, upper_bound, intermediate_results, atol=atol
+                    upper_bound = xp.asarray(self.M if self.closed else self.M - 1, dtype=s.dtype)
+
+            _item_set(
+                results,
+                i,
+                self._length_to_parameter_recursion(
+                    s[i], current_value, lower_bound, upper_bound, xp, intermediate_results, atol=atol
+                ),
+                xp,
             )
 
             lower_bound = results[i]
@@ -1284,9 +1372,9 @@ def curvilinear_reparametrisation_energy(self, atol=1e-6, rtol=1e-6, c=None):
 
         .. _scipy.integrate.quad: https://docs.scipy.org/doc/scipy-1.14.0/reference/generated/scipy.integrate.quad.html
         """
-        arc_length = self.arc_length()
         scale_invariant = c is None
         if scale_invariant:
+            arc_length = self.arc_length()
             c = (arc_length / self.M) ** 2
         upper_limit = self.M if self.closed else self.M - 1
         integral = scipy.integrate.quad(
@@ -1339,9 +1427,9 @@ def curvature(self, t):
         >>> spline = splinebox.Spline(M=M, basis_function=splinebox.Exponential(M), closed=False)
         >>> spline.knots = np.array([[0, 0], [1, 1], [2, 0], [2.5, -0.5], [3, 0]])
         >>> spline.curvature([1, 3])
-        array([ 2.234, -4.104])
+        array([ 2.035, -3.616])
         """
-        self._check_control_points()
+        self._get_control_points_xp_and_device()
         t, single_value = self._convert_to_array(t)
         first_deriv = self(t, derivative=1)
         second_deriv = self(t, derivative=2)
@@ -1414,8 +1502,8 @@ def normal(self, t, frame="bishop", initial_vector=None, shift=None):
         >>> plt.arrow(point[0], point[1], normal[0], normal[1])  # doctest: +SKIP
         >>> plt.show()  # doctest: +SKIP
         """
-        self._check_control_points()
-        t, single_value = self._convert_to_array(t)
+        xp, device = self._get_control_points_xp_and_device()
+        t, single_value = self._convert_to_array(t, xp, device)
         if self.ndim not in (2, 3):
             raise NotImplementedError(
                 f"The normal vector is only implemented for curves in 2D and 3D. Your spline's codomain is {self.ndim} dimensional."
@@ -1537,68 +1625,101 @@ def moving_frame(self, t, method="frenet", initial_vector=None, shift=None):
                 [-1.,  0.,  0.],
                 [ 0., -0.,  1.]]])
         """
-        self._check_control_points()
+        xp, device = self._get_control_points_xp_and_device()
 
         if self.ndim != 3:
             raise RuntimeError("A frame can only be computed for splines in 3D.")
 
-        t, single_value = self._convert_to_array(t)
+        t, single_value = self._convert_to_array(t, xp, device)
 
         # Sort t and keep the indices so the original order can be restored
-        sort_indices = np.argsort(t)
+        sort_indices = xp.argsort(t)
         t = t[sort_indices]
 
         first_derivative = self(t, derivative=1)
 
-        nan_mask = np.isnan(first_derivative)
-        if np.any(nan_mask) and shift is not None:
-            t_mask = np.any(nan_mask, axis=-1)
+        nan_mask = xp.isnan(first_derivative)
+        if xp.any(nan_mask) and shift is not None:
+            t_mask = xp.any(nan_mask, axis=-1)
             first_derivative[t_mask] = self(t[t_mask] + shift, derivative=1)
-        nan_mask = np.isnan(first_derivative)
-        if np.any(nan_mask):
+        nan_mask = xp.isnan(first_derivative)
+        if xp.any(nan_mask):
             raise RuntimeError(
-                f"The frame cannot be computed for t={t[np.any(nan_mask, axis=-1)]} because the spline is not differentiable at those positions. Consider using the `shift` to avoid this problem."
+                f"The frame cannot be computed for t={t[xp.any(nan_mask, axis=-1)]} because the spline is not differentiable at those positions. Consider using the `shift` to avoid this problem."
             )
 
-        frame = np.zeros((len(t), 3, 3))
-        frame[:, 0] = first_derivative / np.linalg.norm(first_derivative, axis=-1)[:, np.newaxis]
+        dtype = (
+            xp.asarray(1.0).dtype
+            if xp.isdtype(self.control_points.dtype, kind="integral") and xp.isdtype(t.dtype, kind="integral")
+            else xp.result_type(self.control_points.dtype, t.dtype)
+        )
+
+        atol = 1e-4 if dtype == xp.float32 else 1e-8
+
+        # Zero value for comparison. pytorch requires a tensor for comparison.
+        zero = xp.asarray(0.0, dtype=dtype, device=device)
+
+        frame = xp.zeros((len(t), 3, 3), dtype=dtype, device=device)
+        frame = _item_set(
+            frame,
+            (slice(len(t)), 0),
+            first_derivative / xp.linalg.vector_norm(first_derivative, axis=-1)[:, xp.newaxis],
+            xp,
+        )
 
         if method == "frenet":
             second_derivative = self(t, derivative=2)
 
-            nan_mask = np.isnan(second_derivative)
-            if np.any(nan_mask) and shift is not None:
-                t_mask = np.any(nan_mask, axis=-1)
+            nan_mask = xp.isnan(second_derivative)
+            if xp.any(nan_mask) and shift is not None:
+                t_mask = xp.any(nan_mask, axis=-1)
                 second_derivative[t_mask] = self(t[t_mask] + shift, derivative=2)
-            nan_mask = np.isnan(second_derivative)
-            if np.any(nan_mask):
+            nan_mask = xp.isnan(second_derivative)
+            if xp.any(nan_mask):
                 raise RuntimeError(
-                    f"The Frenet frame cannot be computed for t={t[np.any(nan_mask, axis=-1)]} because the spline is not twice differentiable at those positions. Consider using `shift` or the Bishop frame."
+                    f"The Frenet frame cannot be computed for t={t[xp.any(nan_mask, axis=-1)]} because the spline is not twice differentiable at those positions. Consider using `shift` or the Bishop frame."
                 )
 
-            frame[:, 2] = np.cross(first_derivative, second_derivative)
-            norm_binormal = np.linalg.norm(frame[:, 2], axis=-1)[:, np.newaxis]
+            frame = _item_set(
+                frame,
+                (slice(len(t)), 2),
+                xp.cross(first_derivative, second_derivative),
+                xp,
+            )
+            norm_binormal = xp.linalg.vector_norm(frame[:, 2], axis=-1)[:, xp.newaxis]
 
-            if np.isclose(norm_binormal[0], 0) or np.isclose(norm_binormal[-1], 0):
+            if xp.isclose(norm_binormal[0], zero, atol=atol) or xp.isclose(norm_binormal[-1], zero, atol=atol):
                 raise RuntimeError(
                     "The Frenet frame cannot be computed at one or both ends of the spline. This is often due to edge padding of the knots. Try to skip t=0 and t=M-1 or change the padding."
                 )
-            if np.any(np.isclose(norm_binormal, 0)):
+            if xp.any(xp.isclose(norm_binormal, zero, atol=atol)):
                 raise RuntimeError(
                     "The Frenet frame is not defined for splines with inflection points or straight segments, try the Bishop frame instead."
                 )
 
-            frame[:, 2] /= norm_binormal
-            frame[:, 1] = np.cross(frame[:, 2], frame[:, 0])
+            frame = _item_set(
+                frame,
+                (slice(len(t)), 2),
+                frame[:, 2] / norm_binormal,
+                xp,
+            )
+            frame = _item_set(
+                frame,
+                (slice(len(t)), 1),
+                xp.cross(frame[:, 2], frame[:, 0]),
+                xp,
+            )
 
         elif method == "bishop":
             if initial_vector is None:
                 tangent = frame[0, 0]
                 # Try to do the same as for the Frenet frame
-                initial_vector = np.cross(np.cross(tangent, self(t[0], derivative=2)), tangent)
-                if np.isclose(np.linalg.norm(initial_vector), 0) or np.any(np.isnan(initial_vector)):
-                    initial_vector = np.zeros(3)
-                    max_axis = np.argmax(np.abs(tangent))
+                initial_vector = xp.cross(xp.cross(tangent, self(t[0].item(), derivative=2)), tangent)
+                if xp.isclose(xp.linalg.vector_norm(initial_vector), zero, atol=atol) or xp.any(
+                    xp.isnan(initial_vector)
+                ):
+                    initial_vector = xp.zeros(3, device=device)
+                    max_axis = xp.argmax(xp.abs(tangent))
                     other_axis = (max_axis + 1) % 3
                     initial_vector[max_axis] = tangent[other_axis]
                     initial_vector[other_axis] = -tangent[max_axis]
@@ -1608,35 +1729,48 @@ def moving_frame(self, t, method="frenet", initial_vector=None, shift=None):
                     UserWarning,
                     stacklevel=2,
                 )
-            initial_vector /= np.linalg.norm(initial_vector)
-            if not np.isclose(np.dot(frame[0, 0], initial_vector), 0):
+            initial_vector = xp.astype(initial_vector, dtype)
+            initial_vector = initial_vector / xp.linalg.vector_norm(initial_vector)
+            if not xp.isclose(xp.dot(frame[0, 0], initial_vector), zero, atol=atol):
                 raise ValueError("The initial vector has to be orthogonal to the tangent at t[0].")
-            frame[0, 1] = initial_vector
-            frame[0, 2] = np.cross(frame[0, 0], initial_vector)
+            frame = _item_set(
+                frame,
+                (0, 1),
+                initial_vector,
+                xp,
+            )
+            frame = _item_set(
+                frame,
+                (0, 2),
+                xp.cross(frame[0, 0], initial_vector),
+                xp,
+            )
             for i in range(1, len(t)):
-                n = np.cross(frame[i - 1, 0], frame[i, 0])
-                norm_n = np.linalg.norm(n)
-                if np.isclose(norm_n, 0):
-                    frame[i, 1] = frame[i - 1, 1]
-                    frame[i, 2] = frame[i - 1, 2]
+                n = xp.cross(frame[i - 1, 0], frame[i, 0])
+                norm_n = xp.linalg.vector_norm(n)
+                if xp.isclose(norm_n, zero, atol=atol):
+                    frame = _item_set(frame, (i, 1), frame[i - 1, 1], xp)
+                    frame = _item_set(frame, (i, 2), frame[i - 1, 2], xp)
                 else:
                     n /= norm_n
-                    phi = np.arccos(np.dot(frame[i - 1, 0], frame[i, 0]))
-                    frame[i, 1] = (
-                        frame[i - 1, 1] * np.cos(phi)
-                        + np.cross(n, frame[i - 1, 1]) * np.sin(phi)
-                        + n * np.dot(n, frame[i - 1, 1]) * (1 - np.cos(phi))
+                    phi = xp.arccos(xp.dot(frame[i - 1, 0], frame[i, 0]))
+                    v1 = (
+                        frame[i - 1, 1] * xp.cos(phi)
+                        + xp.cross(n, frame[i - 1, 1]) * xp.sin(phi)
+                        + n * xp.dot(n, frame[i - 1, 1]) * (1 - xp.cos(phi))
                     )
-                    frame[i, 2] = (
-                        frame[i - 1, 2] * np.cos(phi)
-                        + np.cross(n, frame[i - 1, 2]) * np.sin(phi)
-                        + n * np.dot(n, frame[i - 1, 2]) * (1 - np.cos(phi))
+                    v2 = (
+                        frame[i - 1, 2] * xp.cos(phi)
+                        + xp.cross(n, frame[i - 1, 2]) * xp.sin(phi)
+                        + n * xp.dot(n, frame[i - 1, 2]) * (1 - xp.cos(phi))
                     )
+                    frame = _item_set(frame, (i, 1), v1, xp)
+                    frame = _item_set(frame, (i, 2), v2, xp)
         else:
             raise ValueError(f"Unknown method '{method}' for moving frame.")
 
         # Restore to the original order of t
-        frame = frame[0] if single_value else frame[np.argsort(sort_indices)]
+        frame = frame[0] if single_value else frame[xp.argsort(sort_indices)]
         return frame
 
     @staticmethod
@@ -1656,7 +1790,7 @@ def moving_frame(self, t, method="frenet", initial_vector=None, shift=None):
         nopython=True,
         cache=True,
     )
-    def _compute_tval_and_indices(t, shift, closed, M, pad, tval, indices):
+    def _compute_tval_and_indices_numba(t, shift, closed, M, pad, tval, indices):
         """
         Helper function for creating a csr sparse matrix.
 
@@ -1686,6 +1820,42 @@ def _compute_tval_and_indices(t, shift, closed, M, pad, tval, indices):
         else:
             indices[:] = (t - t_mod_1)[:, np.newaxis] + shift[np.newaxis, :] + pad
 
+    @staticmethod
+    def _compute_tval_and_indices_xp(t, shift, closed, M, pad, xp):
+        """
+        Helper function for computing the t values and indices for
+        non-numpy array backends, mirroring
+        :meth:`splinebox.spline_curves.Spline._compute_tval_and_indices_numba`.
+
+        Parameters
+        ----------
+        t : array
+            The t values at which the spline should be evaluated.
+        shift : np.array
+            The integer displacements that for the control points that affect
+            a given t value. This depends on the half support of the basis function.
+        closed : boolean
+            Whether the spline is closed or not.
+        M : int
+            The number of knots.
+        pad : int
+            The amount of padding the spline has.
+        xp : module
+            The array namespace of the control points.
+
+        Returns
+        -------
+        tval : array
+            The t values at which the basis functions have to be evaluated.
+        indices : array
+            The indices of the control points contributing to each t value.
+        """
+        t_mod_1 = t % 1
+        tval = t_mod_1[:, xp.newaxis] - shift[xp.newaxis, :]
+        indices = xp.asarray(t - t_mod_1, dtype=int)[:, xp.newaxis] + shift[xp.newaxis, :]
+        indices = indices % M if closed else indices + pad
+        return tval, indices
+
     def basis_matrix(self, t, derivative=0):
         r"""
         Computes the basis matrix :math:`\mathbf{\Phi}` as defined in :ref:`theory/data_approximation:Data approximation`.
@@ -1701,51 +1871,77 @@ def basis_matrix(self, t, derivative=0):
         Returns
         -------
         bm : scipy.sparse.csr_array
-            The basis matrix.
+            The basis matrix. If the control points are not numpy arrays, a
+            dense array in the same namespace as the control points is
+            returned instead, since scipy.sparse only supports numpy.
         """
-        t, single_value = self._convert_to_array(t)
+        xp, device = self._get_control_points_xp_and_device()
+        t, single_value = self._convert_to_array(t, xp, device)
 
-        minimum = np.min(t)
-        maximum = np.max(t)
-        if np.isnan(maximum):
+        dtype = (
+            xp.asarray(1.0).dtype
+            if xp.isdtype(t.dtype, kind="integral") and xp.isdtype(self.control_points.dtype, kind="integral")
+            else xp.result_type(t.dtype, self.control_points.dtype)
+        )
+
+        minimum = xp.min(t)
+        maximum = xp.max(t)
+        if bool(xp.isnan(maximum)):
             raise ValueError("t should not contain any NaN values.")
 
         bound = math.ceil(self.half_support)
-        shift = np.arange(-bound + 1, bound + 1)
+        shift = xp.arange(-bound + 1, bound + 1, device=device)
         n_control_points = self.M if self.closed else self.M + 2 * self.pad
 
-        tval = np.empty((len(t), len(shift)), dtype=float)
-        indices = np.empty((len(t), len(shift)), dtype=int)
-        self._compute_tval_and_indices(t, shift, self.closed, self.M, self.pad, tval, indices)
-        data = self.basis_function(tval, derivative=derivative)
+        if array_api_compat.is_numpy_namespace(xp):
+            tval = xp.empty((len(t), len(shift)), dtype=dtype, device=device)
+            indices = xp.empty((len(t), len(shift)), dtype=int, device=device)
+            self._compute_tval_and_indices_numba(t, shift, self.closed, self.M, self.pad, tval, indices)
 
-        if not self.closed and (minimum < 0 or maximum > self.M - 1):
-            mask = (indices >= 0) & (indices <= (self.M - 1 + 2 * self.pad))
+            data = self.basis_function(tval, derivative=derivative)
 
-            # Mask rows that are completely empty
-            row_mask = np.any(mask, axis=1)
-            data = data[row_mask]
-            indices = indices[row_mask]
+            if not self.closed and (minimum < 0 or maximum > self.M - 1):
+                mask = (indices >= 0) & (indices <= (self.M - 1 + 2 * self.pad))
 
-            mask = mask[row_mask]
-            data[~mask] = 0
+                # Mask rows that are completely empty
+                row_mask = np.any(mask, axis=1)
+                data = data[row_mask]
+                indices = indices[row_mask]
 
-            data = data.flatten()
-            indices = indices.flatten()
+                mask = mask[row_mask]
+                data[~mask] = 0
+
+                data = data.flatten()
+                indices = indices.flatten()
 
-            indptr = np.full(len(t) + 1, len(indices))
-            indptr[: len(indices) // len(shift) + 1] = np.arange(0, len(indices) + 1, len(shift))
+                indptr = np.full(len(t) + 1, len(indices))
+                indptr[: len(indices) // len(shift) + 1] = np.arange(0, len(indices) + 1, len(shift))
 
+            else:
+                data = data.flatten()
+                indices = indices.flatten()
+                indptr = np.arange(0, len(indices) + 1, len(shift))
+
+            # This is necessary because there is not enough padding of the control points.
+            # When t=M-1 shift reaches past the last padded control point.
+            # The seems to be the limiting factor. An alternative solution should be found.
+            bm = scipy.sparse.csr_array((data, indices, indptr), shape=(len(t), n_control_points + 1), copy=False)
+            bm = bm[:, :n_control_points]
         else:
-            data = data.flatten()
+            tval, indices = self._compute_tval_and_indices_xp(t, shift, self.closed, self.M, self.pad, xp)
+            tval = xp.astype(tval, dtype)
+            t_indices = xp.tile(xp.arange(len(t))[:, xp.newaxis], (1, tval.shape[1]))
+            tval = tval.flatten()
             indices = indices.flatten()
-            indptr = np.arange(0, len(indices) + 1, len(shift))
+            t_indices = t_indices.flatten()
+            mask = (indices >= 0) & (indices <= (self.M - 1 + 2 * self.pad))
+            tval = tval[mask]
+            indices = indices[mask]
+            t_indices = t_indices[mask]
+            data = self.basis_function(tval, derivative=derivative)
+            bm = xp.zeros((len(t), n_control_points), dtype=data.dtype, device=device)
+            bm = _item_set(bm, (t_indices, indices), data, xp)
 
-        # This is necessary because there is not enough padding of the control points.
-        # When t=M-1 shift reaches past the last padded control point.
-        # The seems to be the limiting factor. An alternative solution should be found.
-        bm = scipy.sparse.csr_array((data, indices, indptr), shape=(len(t), n_control_points + 1), copy=False)
-        bm = bm[:, :n_control_points]
         return bm
 
     def __call__(self, t, derivative=0):
@@ -1780,12 +1976,18 @@ def __call__(self, t, derivative=0):
                [ 1.5,  0.5],
                [ 3. ,  1. ]])
         """
-        self._check_control_points()
-        t, single_value = self._convert_to_array(t)
+        xp, device = self._get_control_points_xp_and_device()
+        t, single_value = self._convert_to_array(t, xp, device)
 
         bm = self.basis_matrix(t, derivative=derivative)
 
-        values = bm @ self.control_points
+        if array_api_compat.is_torch_namespace(xp) and bm.dtype != self.control_points.dtype:
+            # Torch doesn't allow matrix multiplication with mixed dtypes
+            # So we do explicit type-promotion
+            dtype = xp.result_type(bm.dtype, self.control_points.dtype)
+            values = xp.astype(bm, dtype) @ xp.astype(self.control_points, dtype)
+        else:
+            values = bm @ self.control_points
 
         if single_value:
             values = values[0]
@@ -1891,7 +2093,7 @@ def derivative_of_curvature_energy_wrt_control_points(self):
         """
         return 2 * np.einsum("mx,lm->lx", self.control_points, self._h3)
 
-    def _convert_to_array(self, t):
+    def _convert_to_array(self, t, xp, device):
         """
         Helper function that converts a function input
         to an array. This allows functions to accept int and float
@@ -1901,14 +2103,25 @@ def _convert_to_array(self, t):
         final result should be converted back into a single value.
         """
         single_value = False
+        dtype = (
+            self.control_points.dtype
+            if self.control_points is not None and xp.isdtype(self.control_points.dtype, kind="real floating")
+            else xp.asarray(1.0).dtype
+        )
         if not isinstance(t, collections.abc.Iterable):
-            t = np.array([t])
+            t = xp.asarray([t], device=device, dtype=dtype)
             single_value = True
-        elif isinstance(t, np.ndarray) and t.shape == ():
+        elif not array_api_compat.is_array_api_obj(t):
+            t = xp.asarray(t, device=device, dtype=dtype)
+        elif array_api_compat.array_namespace(t) is not xp:
+            raise TypeError(
+                f"The control points are a {xp.__name__.split('.')[-1]} array but "
+                "the parameter values you provided are a "
+                f"{array_api_compat.array_namespace(t).__name__.split('.')[-1]} array."
+            )
+        elif t.shape == ():
             # Array with only one element, e.g. np.array(0.)
-            t = np.array([t.item()])
-        elif not isinstance(t, np.ndarray):
-            t = np.array(t)
+            t = xp.asarray([t.item()], device=device, dtype=t.dtype)
         if t.ndim > 1:
             raise ValueError("The parameter array has to be 1D.")
         return t, single_value
@@ -1919,8 +2132,10 @@ def _control_points_centroid(self):
         and :meth:`splinebox.spline_curves.Spline.rotate`.
         Computes the centroid of the coefficients.
         """
-        self._check_control_points()
-        return np.mean(self.control_points, axis=0)
+        xp, device = self._get_control_points_xp_and_device()
+        return xp.mean(
+            self.control_points, axis=0, dtype=xp.result_type(self.control_points.dtype, xp.asarray(1.0).dtype)
+        )
 
     def translate(self, vector):
         """
@@ -1957,7 +2172,7 @@ def translate(self, vector):
         >>> plt.legend()  # doctest: +SKIP
         >>> plt.show()  # doctest: +SKIP
         """
-        self._check_control_points()
+        self._get_control_points_xp_and_device()
         self.control_points = self.control_points + vector
 
     def scale(self, scaling_factor):
@@ -1995,10 +2210,10 @@ def scale(self, scaling_factor):
         >>> plt.legend()  # doctest: +SKIP
         >>> plt.show()  # doctest: +SKIP
         """
-        self._check_control_points()
+        self._get_control_points_xp_and_device()
         centroid = self._control_points_centroid()
         self.translate(-centroid)
-        self.control_points *= scaling_factor
+        self.control_points = self.control_points * scaling_factor
         self.translate(centroid)
 
     def rotate(self, rotation_matrix, centred=True):
@@ -2038,7 +2253,7 @@ def rotate(self, rotation_matrix, centred=True):
         >>> plt.legend()  # doctest: +SKIP
         >>> plt.show()  # doctest: +SKIP
         """
-        self._check_control_points()
+        self._get_control_points_xp_and_device()
         if self.ndim == 1:
             raise RuntimeError("1D splines can not be rotated.")
 
@@ -2104,7 +2319,11 @@ def distance(self, points, return_t=False):
         >>> plt.gca().set_aspect("equal")  # doctest: +SKIP
         >>> plt.show()  # doctest: +SKIP
         """
-        self._check_control_points()
+        xp, device = self._get_control_points_xp_and_device()
+
+        if not array_api_compat.is_numpy_namespace(xp):
+            raise NotImplementedError("The distance method is only implemented for numpy based splines.")
+
         if self.ndim == 1:
             raise RuntimeError("Cannot compute distance for 1D splines.")
 
@@ -2114,27 +2333,27 @@ def distance(self, points, return_t=False):
             single_point = True
 
         max_t = self.M if self.closed else self.M - 1
-        t = np.linspace(0, max_t, self.M * 10)
+        t = xp.linspace(0, max_t, self.M * 10)
         points_on_spline = self(t)
-        distances = np.linalg.norm(points_on_spline[:, np.newaxis] - points[np.newaxis], axis=-1)
-        t_initial = t[np.argmin(distances, axis=0)]
+        distances = xp.linalg.vector_norm(points_on_spline[:, np.newaxis] - points[np.newaxis], axis=-1)
+        t_initial = t[xp.argmin(distances, axis=0)]
 
         def _distance(t, point):
-            return np.linalg.norm(self(t) - point)
+            return xp.linalg.vector_norm(self(t) - point)
 
         def _jac_distance(t, point):
             spline_point = self(t)
             tangent = self(t, derivative=1)
-            tangent = np.nan_to_num(tangent)
-            norm = np.linalg.norm(spline_point - point)
-            dot_product = np.sum((spline_point - point) * tangent)
+            tangent = xp.nan_to_num(tangent)
+            norm = xp.linalg.vector_norm(spline_point - point)
+            dot_product = xp.sum((spline_point - point) * tangent)
             return 2 / norm * dot_product
 
-        min_distances = np.zeros(points.shape[0])
-        min_t = np.zeros(points.shape[0])
+        min_distances = xp.zeros(points.shape[0], device=device)
+        min_t = xp.zeros(points.shape[0], device=device)
         for i, point in enumerate(points):
             result = scipy.optimize.minimize(_distance, t_initial[i], (point,), jac=_jac_distance, bounds=((0, max_t),))
-            min_distances[i] = np.linalg.norm(self(result.x) - point)
+            min_distances[i] = xp.linalg.vector_norm(self(result.x) - point)
             min_t[i] = result.x[0]
 
         if single_point:
@@ -2302,7 +2521,9 @@ def mesh(
         >>> connectivity.shape
         (6480, 4)
         """
-        self._check_control_points()
+        xp, device = self._get_control_points_xp_and_device()
+        if not array_api_compat.is_numpy_namespace(xp):
+            raise NotImplementedError("Meshes are only implemented for numpy.")
         if self.ndim != 3:
             raise NotImplementedError("Meshes are only implemented for splines in 3D.")
         cap_ends = self._normalize_cap_ends(cap_ends) if mesh_type == "surface" and not self.closed else None
@@ -2561,7 +2782,6 @@ class HermiteSpline(Spline):
     """
 
     _coef_tangent_mismatch_msg = "It looks like control_points and tangents have different shapes."
-    _no_tangents_msg = "This spline object doesn't have any tangents yet."
 
     def __init__(
         self, M, basis_function, closed=False, control_points=None, tangents=None, padding_function=padding_function
@@ -2569,21 +2789,65 @@ def __init__(
         super().__init__(M, basis_function, closed, control_points=control_points, padding_function=padding_function)
         self.tangents = tangents
 
-    def _check_control_points_and_tangents(self):
+    def _get_control_points_and_tangents_xp_and_device(self):
         """
         Most methods require control points and tangents to be set before they
         can be used. This helper function checks if control points and tangents have been
-        set.
+        set and returns their array namespace and the device.
         """
-        self._check_control_points()
+        xp, device = self._get_control_points_xp_and_device()
         if self.tangents is None:
-            raise RuntimeError(self._no_tangents_msg)
+            raise RuntimeError("This spline object doesn't have any tangents yet.")
+        if array_api_compat.array_namespace(self.tangents) is not xp:
+            raise TypeError(
+                f"The control points are a {xp.__name__.split('.')[-1]} array but the tangents are a "
+                f"{array_api_compat.array_namespace(self.tangents).__name__.split('.')[-1]} array."
+            )
+        if array_api_compat.device(self.tangents) != device:
+            raise TypeError(
+                f"Control points are on device {device} and tangents are on device {array_api_compat.device(self.tangents)}. "
+                "They have to be on the same device."
+            )
+        return xp, device
 
     def __repr__(self):
-        return f"splinebox.spline_curves.HermiteSpline(M={repr(self.M)}, basis_function={repr(self.basis_function)}, closed={repr(self.closed)}, control_points=np.{repr(self.control_points)}, tangents=np.{repr(self.tangents)})"
+        xp, device = self._get_control_points_xp_and_device()
+
+        control_points_repr = repr(self.control_points)
+        tangents_repr = repr(self.tangents)
+
+        control_points_repr = (
+            re.split(r",\s*dtype\b", control_points_repr)[0]
+            if "dtype" in control_points_repr
+            else control_points_repr[:-1]
+        )
+        tangents_repr = re.split(r",\s*dtype\b", tangents_repr)[0] if "dtype" in tangents_repr else tangents_repr[:-1]
+
+        dtype_repr = f"{xp.__name__}.{repr(self.control_points.dtype)}"
+
+        if array_api_compat.is_numpy_namespace(xp):
+            control_points_repr = f"np.{control_points_repr}, dtype={dtype_repr})"
+            tangents_repr = f"np.{tangents_repr}, dtype={dtype_repr})"
+        elif array_api_compat.is_jax_namespace(xp):
+            device_repr = f"jax.devices('{device.platform}')[{device.id}]"
+            control_points_repr = control_points_repr.replace("Array", "array")
+            control_points_repr = f"{xp.__name__}.{control_points_repr}, dtype={dtype_repr}, device={device_repr})"
+            tangents_repr = tangents_repr.replace("Array", "array")
+            tangents_repr = f"{xp.__name__}.{tangents_repr}, dtype={dtype_repr}, device={device_repr})"
+        else:
+            device_repr = f"{xp.__name__}.{repr(device)}"
+            control_points_repr = f"{xp.__name__}.{control_points_repr}, dtype={dtype_repr}, device={device_repr})"
+            tangents_repr = f"{xp.__name__}.{tangents_repr}, dtype={dtype_repr}, device={device_repr})"
+
+        return f"splinebox.spline_curves.HermiteSpline(M={repr(self.M)}, basis_function={repr(self.basis_function)}, closed={repr(self.closed)}, control_points={control_points_repr}, tangents={tangents_repr})"
 
     def __eq__(self, other):
-        return super().__eq__(other) and np.all(self.tangents == other.tangents)
+        xp, device = self._get_control_points_and_tangents_xp_and_device()
+        return (
+            super().__eq__(other) and xp.all(self.tangents == other.tangents)
+            if self.tangents is not None and other.tangents is not None
+            else self.tangents is other.tangents
+        )
 
     @property
     def tangents(self):
@@ -2624,16 +2888,23 @@ def basis_function(self, value):
         self._basis_function = value
 
     def fit(self, points, boundary_condition="free"):
-        if len(points) < 2 * (self.M + 2 * self.pad):
+        if not array_api_compat.is_array_api_obj(points):
+            raise TypeError("points has to be an array api compatible object, e.g. numpy array, torch tensor, ...")
+        xp = array_api_compat.array_namespace(points)
+        device = array_api_compat.device(points)
+
+        n_points = len(points)
+        if n_points < 2 * (self.M + 2 * self.pad):
             raise RuntimeError(
                 f"You provided too few points. For a unique solution you need to provide at least 2*({self.M}+{2 * self.pad}) points to match the number of control points and tangents (including padding). You provided {len(points)} points. Consider providing more points or reducing the number of knots M."
             )
-        if len(points) < self.M:
-            raise RuntimeError(
-                "You provided fewer data points than your spline has knots. For the fit to have a unique solution you need to provide at least as many data points as your spline has knots. Consider adding more data or reducing the number of knots M."
+        if points.ndim != 2:
+            raise ValueError(
+                "The points should be provided as a 2D array, with the first dimension encoding the different data points and the second dimension encoding the dimensionality of the data points."
             )
 
-        n_points = len(points)
+        # We do the fit in numpy to be able to exploit scipy's sparse arrays
+        points = np.array(points)
         n_control_points = self.M if self.closed else self.M + 2 * self.pad
 
         t = np.linspace(0, self.M, n_points + 1)[:-1] if self.closed else np.linspace(0, self.M - 1, n_points)
@@ -2642,7 +2913,7 @@ def fit(self, points, boundary_condition="free"):
         shift = np.arange(-bound + 1, bound + 1)
         tval = np.empty((len(t), len(shift)), dtype=float)
         indices = np.empty((len(t), len(shift)), dtype=int)
-        self._compute_tval_and_indices(t, shift, self.closed, self.M, self.pad, tval, indices)
+        self._compute_tval_and_indices_numba(t, shift, self.closed, self.M, self.pad, tval, indices)
 
         row = np.repeat(np.arange(indices.shape[0]), indices.shape[1] * 2)
         col = np.empty(2 * indices.shape[0] * indices.shape[1])
@@ -2666,10 +2937,8 @@ def fit(self, points, boundary_condition="free"):
                 self.control_points = solution[:half]
                 self.tangents = solution[half:]
             else:
-                if self.control_points is None:
-                    self.control_points = np.empty((n_control_points, points.shape[1]))
-                if self.tangents is None:
-                    self.tangents = np.empty((n_control_points, points.shape[1]))
+                self.control_points = np.empty((n_control_points, points.shape[1]))
+                self.tangents = np.empty((n_control_points, points.shape[1]))
                 for i in range(self.ndim):
                     solution = scipy.sparse.linalg.lsqr(basis_function_values, points[:, i])[0]
                     self.control_points[:, i] = solution[:half]
@@ -2693,12 +2962,10 @@ def fit(self, points, boundary_condition="free"):
                 (data, (row, col)), shape=(n_points, 2 * n_control_points - 2)
             )
 
-            if self.control_points is None:
-                if points.ndim == 1:
-                    points = points[:, np.newaxis]
-                self.control_points = np.empty((n_control_points, points.shape[1]))
-            if self.tangents is None:
-                self.tangents = np.empty((n_control_points, points.shape[1]))
+            if points.ndim == 1:
+                points = points[:, np.newaxis]
+            self.control_points = np.empty((n_control_points, points.shape[1]))
+            self.tangents = np.empty((n_control_points, points.shape[1]))
             for i in range(self.ndim):
                 solution = scipy.sparse.linalg.lsqr(basis_function_values, points[:, i])[0]
                 self.control_points[:, i] = solution[:n_control_points]
@@ -2715,64 +2982,104 @@ def fit(self, points, boundary_condition="free"):
         else:
             raise ValueError(f"Unknown boundary_conditions {boundary_condition}")
 
+        # Convert to the correct array type
+        if not array_api_compat.is_numpy_namespace(xp):
+            self.control_points = xp.asarray(self.control_points, device=device)
+            self.tangents = xp.asarray(self.tangents, device=device)
+
     def basis_matrix(self, t, derivative=0):
-        t, single_value = self._convert_to_array(t)
+        xp, device = self._get_control_points_and_tangents_xp_and_device()
+        t, single_value = self._convert_to_array(t, xp, device)
+
+        dtype = (
+            xp.asarray(1.0).dtype
+            if xp.isdtype(t.dtype, kind="integral") and xp.isdtype(self.control_points.dtype, kind="integral")
+            else xp.result_type(t.dtype, self.control_points.dtype)
+        )
 
-        minimum = np.min(t)
-        maximum = np.max(t)
-        if np.isnan(maximum):
+        minimum = xp.min(t)
+        maximum = xp.max(t)
+        if bool(xp.isnan(maximum)):
             raise ValueError("t should not contain any NaN values.")
 
         bound = math.ceil(self.half_support)
-        shift = np.arange(-bound + 1, bound + 1)
+        shift = xp.arange(-bound + 1, bound + 1, device=device)
         n_control_points = self.M if self.closed else self.M + 2 * self.pad
 
-        tval = np.empty((len(t), len(shift)), dtype=float)
-        indices = np.empty((len(t), len(shift)), dtype=int)
-        self._compute_tval_and_indices(t, shift, self.closed, self.M, self.pad, tval, indices)
-        data = self.basis_function(tval, derivative=derivative)
-        data0 = data[..., 0]
-        data1 = data[..., 1]
+        if array_api_compat.is_numpy_namespace(xp):
+            tval = xp.empty((len(t), len(shift)), dtype=dtype, device=device)
+            indices = xp.empty((len(t), len(shift)), dtype=int, device=device)
+            self._compute_tval_and_indices_numba(t, shift, self.closed, self.M, self.pad, tval, indices)
+            data = self.basis_function(tval, derivative=derivative)
+            data0 = data[..., 0]
+            data1 = data[..., 1]
 
-        if not self.closed and (minimum < 0 or maximum > self.M - 1):
-            mask = (indices >= 0) & (indices <= (self.M - 1 + 2 * self.pad))
+            if not self.closed and (minimum < 0 or maximum > self.M - 1):
+                mask = (indices >= 0) & (indices <= (self.M - 1 + 2 * self.pad))
 
-            # Mask rows that are completely empty
-            row_mask = np.any(mask, axis=1)
-            data0 = data0[row_mask]
-            data1 = data1[row_mask]
-            indices = indices[row_mask]
+                # Mask rows that are completely empty
+                row_mask = np.any(mask, axis=1)
+                data0 = data0[row_mask]
+                data1 = data1[row_mask]
+                indices = indices[row_mask]
 
-            mask = mask[row_mask]
-            data0[~mask] = 0
-            data1[~mask] = 0
+                mask = mask[row_mask]
+                data0[~mask] = 0
+                data1[~mask] = 0
 
-            data0 = data0.flatten()
-            data1 = data1.flatten()
-            indices = indices.flatten()
+                data0 = data0.flatten()
+                data1 = data1.flatten()
+                indices = indices.flatten()
 
-            indptr = np.full(len(t) + 1, len(indices))
-            indptr[: len(indices) // len(shift) + 1] = np.arange(0, len(indices) + 1, len(shift))
+                indptr = np.full(len(t) + 1, len(indices))
+                indptr[: len(indices) // len(shift) + 1] = np.arange(0, len(indices) + 1, len(shift))
 
+            else:
+                data0 = data0.flatten()
+                data1 = data1.flatten()
+                indices = indices.flatten()
+                indptr = np.arange(0, len(indices) + 1, len(shift))
+            bm0 = scipy.sparse.csr_array((data0, indices, indptr), shape=(len(t), n_control_points + 1), copy=False)
+            bm0 = bm0[:, :n_control_points]
+            bm1 = scipy.sparse.csr_array((data1, indices, indptr), shape=(len(t), n_control_points + 1), copy=False)
+            bm1 = bm1[:, :n_control_points]
         else:
-            data0 = data0.flatten()
-            data1 = data1.flatten()
+            tval, indices = self._compute_tval_and_indices_xp(t, shift, self.closed, self.M, self.pad, xp)
+            tval = xp.astype(tval, dtype)
+            t_indices = xp.tile(xp.arange(len(t))[:, xp.newaxis], (1, tval.shape[1]))
+            tval = tval.flatten()
             indices = indices.flatten()
-            indptr = np.arange(0, len(indices) + 1, len(shift))
-
-        bm0 = scipy.sparse.csr_array((data0, indices, indptr), shape=(len(t), n_control_points + 1), copy=False)
-        bm0 = bm0[:, :n_control_points]
-        bm1 = scipy.sparse.csr_array((data1, indices, indptr), shape=(len(t), n_control_points + 1), copy=False)
-        bm1 = bm1[:, :n_control_points]
+            t_indices = t_indices.flatten()
+            mask = (indices >= 0) & (indices <= (self.M - 1 + 2 * self.pad))
+            tval = tval[mask]
+            indices = indices[mask]
+            t_indices = t_indices[mask]
+            data = self.basis_function(tval, derivative=derivative)
+            data0 = data[..., 0]
+            data1 = data[..., 1]
+            bm0 = xp.zeros((len(t), n_control_points), dtype=data.dtype, device=device)
+            bm1 = xp.zeros((len(t), n_control_points), dtype=data.dtype, device=device)
+            bm0 = _item_set(bm0, (t_indices, indices), data0, xp)
+            bm1 = _item_set(bm1, (t_indices, indices), data1, xp)
         return bm0, bm1
 
     def __call__(self, t, derivative=0):
-        self._check_control_points_and_tangents()
-        t, single_value = self._convert_to_array(t)
+        xp, device = self._get_control_points_and_tangents_xp_and_device()
+        t, single_value = self._convert_to_array(t, xp, device)
 
         bm0, bm1 = self.basis_matrix(t, derivative=derivative)
 
-        values = bm0 @ self.control_points + bm1 @ self.tangents
+        if array_api_compat.is_torch_namespace(xp) and (
+            bm0.dtype != self.control_points.dtype or bm1.dtype != self.tangents.dtype
+        ):
+            # Torch doesn't allow matrix multiplication with mixed dtypes
+            # So we do explicit type-promotion
+            dtype = xp.result_type(bm0.dtype, self.control_points.dtype, self.tangents.dtype)
+            values = xp.astype(bm0, dtype) @ xp.astype(self.control_points, dtype) + xp.astype(bm1, dtype) @ xp.astype(
+                self.tangents, dtype
+            )
+        else:
+            values = bm0 @ self.control_points + bm1 @ self.tangents
 
         if single_value:
             values = values[0]
@@ -2828,14 +3135,14 @@ def derivative_of_curvature_energy_wrt_control_points(self):
         )
 
     def scale(self, scaling_factor):
-        self._check_control_points_and_tangents()
-        Spline.scale(self, scaling_factor)
-        self.tangents *= scaling_factor
+        xp, device = self._get_control_points_and_tangents_xp_and_device()
+        super().scale(scaling_factor)
+        self.tangents = self.tangents * scaling_factor
 
     def rotate(self, rotation_matrix, centred=True):
-        self._check_control_points_and_tangents()
+        xp, device = self._get_control_points_and_tangents_xp_and_device()
 
-        Spline.rotate(self, rotation_matrix, centred=centred)
+        self.rotate(self, rotation_matrix, centred=centred)
 
         for k in range(len(self.tangents)):
             self.tangents[k] = np.matmul(rotation_matrix, self.tangents[k])
@@ -2902,7 +3209,6 @@ def splines_to_json(path, splines, version=1):
     Create two random splines...
 
     >>> spline1 = splinebox.Spline(M=5, basis_function=splinebox.B3(), closed=False)
-    >>> spline1.control_points = np.random.rand(7, 3)
     >>> spline2 = splinebox.Spline(M=6, basis_function=splinebox.CatmullRom(), closed=True)
     >>> spline2.control_points = np.random.rand(6, 2)
 
diff --git a/tests/conftest.py b/tests/conftest.py
index ab5a6d7..b7bff0a 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,13 +1,95 @@
+import collections
+import contextlib
 import functools
 import itertools
 import math
-import unittest.mock
 
-import numpy as np
 import pytest
 import splinebox.basis_functions
 import splinebox.spline_curves
 
+ARRAY_BACKENDS_AND_DEVICES = []
+
+# NumPy
+with contextlib.suppress(ImportError):
+
+    import array_api_compat.numpy as np
+
+    ARRAY_BACKENDS_AND_DEVICES.append(pytest.param((np, "cpu"), id="numpy-cpu"))
+
+# PyTorch
+with contextlib.suppress(ImportError):
+
+    import array_api_compat.torch as torch
+    import torch as _torch
+
+    ARRAY_BACKENDS_AND_DEVICES.append(pytest.param((torch, torch.device("cpu")), id="torch-cpu"))
+
+    if _torch.cuda.is_available():
+        ARRAY_BACKENDS_AND_DEVICES.append(pytest.param((torch, torch.device("cuda", index=0)), id="torch-cuda"))
+
+# JAX
+with contextlib.suppress(ImportError):
+
+    import jax
+    import jax.numpy
+
+    # Enable 64-bit floats so that jax matches the float64 semantics of the
+    # other array backends (jax silently truncates float64 data otherwise).
+    jax.config.update("jax_enable_x64", True)
+
+    ARRAY_BACKENDS_AND_DEVICES.append(pytest.param((jax.numpy, jax.devices("cpu")[0]), id="jax-cpu"))
+
+    if any(device.platform == "gpu" for device in jax.devices()):
+        gpus = jax.devices("gpu")
+        ARRAY_BACKENDS_AND_DEVICES.append(pytest.param((jax.numpy, gpus[0]), id="jax-gpu"))
+
+
+@pytest.fixture(params=ARRAY_BACKENDS_AND_DEVICES)
+def xp_device(request):
+    return request.param
+
+
+@pytest.fixture
+def xp(xp_device):
+    return xp_device[0]
+
+
+@pytest.fixture
+def device(xp_device):
+    return xp_device[1]
+
+
+@pytest.fixture(params=["int", "float32", "float64"])
+def dtype(request, xp):
+    """
+    The dtype of the input array.
+
+    The output is expected to be of dtype float32 or
+    float64 if the input is of dtype float32 or float64 respectively and the
+    default floating point dtype of the array namespace if the input is of an
+    integer dtype.
+    """
+    return int if request.param == "int" else getattr(xp, request.param)
+
+
+@pytest.fixture
+def as_xp_array(xp, device):
+    def _as_xp_array(obj, **kwargs):
+        return xp.asarray(obj, device=device, **kwargs)
+
+    return _as_xp_array
+
+
+@pytest.fixture(
+    params=[
+        # pytest.param(True, id="jit"),
+        pytest.param(False, id="no-jit"),
+    ]
+)
+def jit(request):
+    return request.param
+
 
 @pytest.fixture(
     params=[
@@ -97,12 +179,12 @@ def _is_locally_refinable(basis_function):
     return _is_locally_refinable
 
 
-@pytest.fixture(params=[4, 5, 8])
+@pytest.fixture(params=[4, 5, 8], ids=["M=4", "M=5", "M=8"])
 def M(request):
     return request.param
 
 
-@pytest.fixture(params=[True, False])
+@pytest.fixture(params=[True, False], ids=["closed", "open"])
 def closed(request):
     return request.param
 
@@ -120,11 +202,8 @@ def spline_curve(basis_function, M, closed, padding_function):
         spline = splinebox.spline_curves.HermiteSpline(
             M, basis_function, closed=closed, padding_function=padding_function
         )
-        spline._check_control_points_and_tangents = unittest.mock.MagicMock()
-        spline._check_control_points = unittest.mock.MagicMock()
     else:
         spline = splinebox.spline_curves.Spline(M, basis_function, closed=closed, padding_function=padding_function)
-        spline._check_control_points = unittest.mock.MagicMock()
     return spline
 
 
@@ -173,17 +252,36 @@ def curvature_energy_compatible_spline_curve(curvature_energy_compatible_basis_f
     return splinebox.spline_curves.Spline(M, curvature_energy_compatible_basis_function, closed=False)
 
 
-@pytest.fixture(params=[0, 1, 2])
+@pytest.fixture(params=[0, 1, 2], ids=["derivative=0", "derivative=1", "derivative=2"])
 def derivative(request):
     return request.param
 
 
 @pytest.fixture(params=[1.5, np.array(0.0), np.linspace(0, 10, 1000)])
-def call_positions(request):
+def call_positions(request, as_xp_array):
+    if isinstance(request.param, collections.abc.Iterable):
+        return as_xp_array(request.param)
+    return request.param
+
+
+@pytest.fixture(
+    params=[
+        0.5,
+        1,
+        4,
+        [0.5],
+        [0, 1, 4],
+        [-0.5, 0.0, 0.5, 4.1],
+        [[-0.7, -0.2, 0.3], [-0.6, -0.1, 0.4], [-0.5, 0.0, 0.5], [-0.4, 0.1, 0.6]],
+    ]
+)
+def basis_function_call_positions(request, as_xp_array, dtype):
+    if isinstance(request.param, collections.abc.Iterable):
+        return as_xp_array(request.param, dtype=dtype)
     return request.param
 
 
-@pytest.fixture(params=[1, 2, 3, 5])
+@pytest.fixture(params=[1, 2, 3, 5], ids=["ndim=1", "ndim=2", "ndim=3", "ndim=5"])
 def codomain_dimensionality(request):
     """
     The number of dimensions the spline output
@@ -194,37 +292,43 @@ def codomain_dimensionality(request):
     return request.param
 
 
-def _control_point_gen(M, support, closed, ndim, rng=None):
-    if rng is None:
-        rng = np.random.default_rng(seed=42)
-    points = rng.random((M, ndim)) * 100 if closed else rng.random((M + 2 * (math.ceil(support / 2) - 1), ndim)) * 100
-    return points
+def _control_point_gen(_M, _support, _closed, _ndim, _as_xp_array, _dtype, _rng=None):
+    if _rng is None:
+        _rng = np.random.default_rng(seed=42)
+    points = (
+        _rng.random((_M, _ndim)) * 100
+        if _closed
+        else _rng.random((_M + 2 * (math.ceil(_support / 2) - 1), _ndim)) * 100
+    )
+    return _as_xp_array(points, dtype=_dtype)
 
 
 @pytest.fixture
-def coeff_gen(codomain_dimensionality):
+def coeff_gen(codomain_dimensionality, as_xp_array, dtype):
     rng = np.random.default_rng(seed=1492)
-    return functools.partial(_control_point_gen, ndim=codomain_dimensionality, rng=rng)
+    return functools.partial(
+        _control_point_gen, _ndim=codomain_dimensionality, _as_xp_array=as_xp_array, _dtype=dtype, _rng=rng
+    )
 
 
 @pytest.fixture
-def control_point_gen_2D():
+def control_point_gen_2D(as_xp_array, dtype):
     rng = np.random.default_rng(seed=9052)
-    return functools.partial(_control_point_gen, ndim=2, rng=rng)
+    return functools.partial(_control_point_gen, _ndim=2, _as_xp_array=as_xp_array, _dtype=dtype, _rng=rng)
 
 
 @pytest.fixture
-def control_point_gen_3D():
+def control_point_gen_3D(as_xp_array, dtype):
     rng = np.random.default_rng(seed=7755)
-    return functools.partial(_control_point_gen, ndim=3, rng=rng)
+    return functools.partial(_control_point_gen, _ndim=3, _as_xp_array=as_xp_array, _dtype=dtype, _rng=rng)
 
 
 @pytest.fixture
-def knot_gen(codomain_dimensionality):
+def knot_gen(codomain_dimensionality, as_xp_array, dtype):
     rng = np.random.default_rng(seed=2657)
 
     def _knot_gen(M=100):
-        return rng.random((M, codomain_dimensionality))
+        return as_xp_array(rng.random((M, codomain_dimensionality)), dtype=dtype)
 
     return _knot_gen
 
@@ -259,11 +363,11 @@ def is_hermite_spline():
 
 
 @pytest.fixture
-def translation_vector(codomain_dimensionality):
-    return np.random.rand(codomain_dimensionality) * 10 - 5
+def translation_vector(codomain_dimensionality, as_xp_array, dtype):
+    return as_xp_array(np.random.rand(codomain_dimensionality) * 10 - 5, dtype=dtype)
 
 
-def _initialize_spline_curve(spline, control_point_gen):
+def _initialise_spline_curve(spline, control_point_gen):
     support = spline.basis_function.support
     closed = spline.closed
     M = spline.M
@@ -276,58 +380,63 @@ def _initialize_spline_curve(spline, control_point_gen):
 
 
 @pytest.fixture
-def initialize_spline_curve(is_hermite_spline, coeff_gen):
-    return functools.partial(_initialize_spline_curve, control_point_gen=coeff_gen)
+def initialise_spline_curve(is_hermite_spline, coeff_gen):
+    return functools.partial(_initialise_spline_curve, control_point_gen=coeff_gen)
 
 
 @pytest.fixture
-def initialize_2D_spline_curve(is_hermite_spline, control_point_gen_2D):
-    return functools.partial(_initialize_spline_curve, control_point_gen=control_point_gen_2D)
+def initialise_2D_spline_curve(is_hermite_spline, control_point_gen_2D):
+    return functools.partial(_initialise_spline_curve, control_point_gen=control_point_gen_2D)
 
 
 @pytest.fixture
-def initialize_3D_spline_curve(is_hermite_spline, control_point_gen_3D):
-    return functools.partial(_initialize_spline_curve, control_point_gen=control_point_gen_3D)
+def initialise_3D_spline_curve(is_hermite_spline, control_point_gen_3D):
+    return functools.partial(_initialise_spline_curve, control_point_gen=control_point_gen_3D)
 
 
 @pytest.fixture
-def initialized_spline_curve(spline_curve, initialize_spline_curve):
-    return initialize_spline_curve(spline_curve)
+def initialised_spline_curve(spline_curve, initialise_spline_curve):
+    return initialise_spline_curve(spline_curve)
 
 
 @pytest.fixture
-def initialized_twice_differentiable_spline_curve(twice_differentiable_spline_curve, initialize_spline_curve):
-    return initialize_spline_curve(twice_differentiable_spline_curve)
+def initialised_twice_differentiable_spline_curve(twice_differentiable_spline_curve, initialise_spline_curve):
+    return initialise_spline_curve(twice_differentiable_spline_curve)
 
 
 @pytest.fixture
-def initialized_2D_spline_curve(spline_curve, initialize_2D_spline_curve):
-    return initialize_2D_spline_curve(spline_curve)
+def initialised_2D_spline_curve(spline_curve, initialise_2D_spline_curve):
+    return initialise_2D_spline_curve(spline_curve)
 
 
 @pytest.fixture
-def initialized_3D_spline_curve(spline_curve, initialize_3D_spline_curve):
-    return initialize_3D_spline_curve(spline_curve)
+def initialised_3D_spline_curve(spline_curve, initialise_3D_spline_curve):
+    return initialise_3D_spline_curve(spline_curve)
 
 
 @pytest.fixture
-def initialized_twice_differentiable_2D_spline_curve(twice_differentiable_spline_curve, initialize_2D_spline_curve):
-    return initialize_2D_spline_curve(twice_differentiable_spline_curve)
+def initialised_twice_differentiable_2D_spline_curve(twice_differentiable_spline_curve, initialise_2D_spline_curve):
+    return initialise_2D_spline_curve(twice_differentiable_spline_curve)
 
 
 @pytest.fixture
-def initialized_twice_differentiable_3D_spline_curve(twice_differentiable_spline_curve, initialize_3D_spline_curve):
-    return initialize_3D_spline_curve(twice_differentiable_spline_curve)
+def initialised_twice_differentiable_3D_spline_curve(twice_differentiable_spline_curve, initialise_3D_spline_curve):
+    return initialise_3D_spline_curve(twice_differentiable_spline_curve)
 
 
 @pytest.fixture
-def initialized_non_hermite_spline_curve(non_hermite_spline_curve, initialize_spline_curve):
-    return initialize_spline_curve(non_hermite_spline_curve)
+def initialised_non_hermite_spline_curve(non_hermite_spline_curve, initialise_spline_curve):
+    return initialise_spline_curve(non_hermite_spline_curve)
 
 
 @pytest.fixture
-def initialized_hermite_spline_curve(hermite_spline_curve, initialize_spline_curve):
-    return initialize_spline_curve(hermite_spline_curve)
+def initialised_hermite_spline_curve(hermite_spline_curve, initialise_spline_curve):
+    return initialise_spline_curve(hermite_spline_curve)
+
+
+@pytest.fixture
+def initialised_closed_spline_curve(closed_spline_curve, initialise_spline_curve):
+    return initialise_spline_curve(closed_spline_curve)
 
 
 @pytest.fixture
@@ -356,15 +465,18 @@ def boundary_condition(request):
     return request.param
 
 
-@pytest.fixture(params=[2, 5, 20, 50])
+@pytest.fixture(
+    params=[2, 5, 20, 50],
+    ids=["2 data points", "5 data points", "20 data points", "50 data points"],
+)
 def n_points(request):
     return request.param
 
 
 @pytest.fixture
-def points(codomain_dimensionality, n_points):
+def points(codomain_dimensionality, n_points, as_xp_array):
     rng = np.random.default_rng(seed=5544)
-    return rng.random((n_points, codomain_dimensionality))
+    return as_xp_array(rng.random((n_points, codomain_dimensionality)))
 
 
 @pytest.fixture(params=[False, True])
@@ -427,5 +539,5 @@ def single_value_t(request):
         np.array([1, 2, 3], dtype=int),
     ]
 )
-def array_t(request):
-    return request.param
+def array_t(request, as_xp_array, dtype):
+    return as_xp_array(request.param, dtype=dtype)
diff --git a/tests/test_basis_functions.py b/tests/test_basis_functions.py
index 5225308..98f5f9b 100644
--- a/tests/test_basis_functions.py
+++ b/tests/test_basis_functions.py
@@ -1,29 +1,33 @@
+import collections
 import math
 
+import array_api_compat
 import numpy as np
 import pytest
 import splinebox.basis_functions
 
 
-def test_base_class_call(derivative):
+def test_base_class_call(derivative, xp, as_xp_array):
     basis_function = splinebox.basis_functions.BasisFunction(False, 2)
-    x = np.arange(10)
+    x = as_xp_array(xp.linspace(-4, 4, 10))
+
     with pytest.raises(NotImplementedError):
-        basis_function(x, derivative=derivative)
+        # The pytorch compiler does not work for raised errors.
+        basis_function(x, derivative=derivative, jit=False)
 
 
-def test_base_class_filters_and_refinement_mask():
+def test_base_class_filters_and_refinement_mask(xp, as_xp_array):
     basis_function = splinebox.basis_functions.BasisFunction(False, 2)
-    s = np.arange(10)
+    s = as_xp_array(xp.linspace(0, 4, 10))
     with pytest.raises(NotImplementedError):
-        basis_function.filter_symmetric(s)
+        basis_function.filter_symmetric(s, jit=False)
     with pytest.raises(NotImplementedError):
-        basis_function.filter_periodic(s)
+        basis_function.filter_periodic(s, jit=False)
     with pytest.raises(NotImplementedError):
         basis_function.refinement_mask()
 
 
-def test_filters(basis_function, is_interpolating, knot_gen, request):
+def test_filters(basis_function, is_interpolating, knot_gen, xp, dtype, request):
     if isinstance(basis_function, splinebox.basis_functions.B2):
         # The filter_symmetric and filter_periodic are not implemented for B2
         request.node.add_marker(pytest.mark.xfail)
@@ -32,19 +36,31 @@ def test_filters(basis_function, is_interpolating, knot_gen, request):
     if is_interpolating(basis_function):
         # The filters should not do anything since knot and control points are
         # the same for interpolating splines.
-        assert np.allclose(basis_function.filter_symmetric(s), s)
-        assert np.allclose(basis_function.filter_periodic(s), s)
-    else:
-        # Flipping the knot or the control points should give the same result
-        assert np.allclose(basis_function.filter_symmetric(s), basis_function.filter_symmetric(s[::-1])[::-1])
-        assert np.allclose(basis_function.filter_periodic(s), basis_function.filter_periodic(s[::-1])[::-1])
+        filtered = basis_function.filter_symmetric(s)
+        assert xp.allclose(filtered, xp.astype(s, filtered.dtype))
+        filtered = basis_function.filter_periodic(s)
+        assert xp.allclose(filtered, xp.astype(s, filtered.dtype))
+
+    # Flipping the knot or the control points should give the same result
+    flipped = xp.flip(s, axis=0)
+    atol = 1e-5 if s.dtype is int or s.dtype == xp.float32 else 1e-8
+    assert xp.allclose(
+        basis_function.filter_symmetric(s),
+        xp.flip(basis_function.filter_symmetric(flipped), axis=0),
+        atol=atol,
+    )
+    assert xp.allclose(
+        basis_function.filter_periodic(s),
+        xp.flip(basis_function.filter_periodic(flipped), axis=0),
+        atol=atol,
+    )
 
 
 def test_refinement_mask(basis_function, is_locally_refinable, request):
     if isinstance(basis_function, splinebox.basis_functions.Exponential):
         # The refinement mask of the exponential basis function fails.
         # This might be because the value of M has to be changed after
-        # refinement. Since we don't have an imidiate need for the refinement
+        # refinement. Since we don't have an immediate need for the refinement
         # mask, we xfail the tests for now and will revisit this later.
         request.node.add_marker(pytest.mark.xfail)
     if is_locally_refinable(basis_function):
@@ -70,48 +86,75 @@ def test_refinement_mask(basis_function, is_locally_refinable, request):
             basis_function.refinement_mask()
 
 
-def test_call_output_type(basis_function, derivative, is_hermite_basis_function):
+def test_call_output_type(
+    basis_function, derivative, basis_function_call_positions, dtype, is_hermite_basis_function, as_xp_array, xp
+):
     """
     If a single value is provided, a single value should be returned.
-    Otherwise, a numpy array should be returned.
+    Otherwise, an array should be returned.
+
+    The value returned should be a float. For integer input the returned
+    value should have the default floating point dtype of the array namespace
+    while float32 and float64 input should be returned unchanged.
     """
-    for t in [
-        0.5,
-        np.array([0.5]),
-        np.array([-0.5, 0.0, 0.5]),
-        np.array([[-0.7, -0.2, 0.3], [-0.6, -0.1, 0.4], [-0.5, 0.0, 0.5], [-0.4, 0.1, 0.6]]),
-    ]:
+    t = basis_function_call_positions
+
+    vals = basis_function(t, derivative=derivative)
 
-        vals = basis_function(t, derivative=derivative)
+    if not isinstance(t, collections.abc.Iterable):
 
         if is_hermite_basis_function(basis_function):
-            assert isinstance(vals, np.ndarray)
-            assert vals.shape[-1] == 2
-            if isinstance(t, float):
-                assert vals.ndim == 1
-            else:
-                assert np.all(vals.shape[:-1] == t.shape)
+            assert vals.ndim == 1
+            assert vals.shape[0] == 2
+        else:
+            assert isinstance(vals, float)
 
+    else:
+
+        # Check that the output is an array from the xp namespace
+        assert array_api_compat.is_array_api_obj(vals)
+        assert array_api_compat.array_namespace(vals) is xp
+
+        # Integer input should be cast to the default floating point dtype of
+        # the array namespace while float32 and float64 input keeps its dtype.
+        assert xp.isdtype(vals.dtype, kind="real floating")
+        expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+        assert vals.dtype == expected_dtype
+
+        if is_hermite_basis_function(basis_function):
+            assert np.all(vals.shape[:-1] == t.shape)
         else:
-            assert isinstance(vals, type(t))
-            if isinstance(t, np.ndarray):
-                assert np.all(vals.shape == t.shape)
+            assert np.all(vals.shape == t.shape)
 
 
-def test_filter_output_type(basis_function, request):
+def test_filter_output_type(basis_function, as_xp_array, dtype, xp, request):
 
     if isinstance(basis_function, splinebox.basis_functions.B2):
         # The filter_symmetric and filter_periodic are not implemented for B2
         request.node.add_marker(pytest.mark.xfail)
 
     for s in [np.array([1, 2, 3, 4, 5]), np.array([[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]])]:
+        s = as_xp_array(s, dtype=dtype)
         output = basis_function.filter_periodic(s)
-        assert isinstance(output, np.ndarray)
-        assert np.issubdtype(output.dtype, np.floating)
+
+        assert array_api_compat.is_array_api_obj(output)
+        assert array_api_compat.array_namespace(output) is xp
+        # Integer input should be cast to the default floating point dtype of
+        # the array namespace while float32 and float64 input keeps its dtype.
+        assert xp.isdtype(output.dtype, kind="real floating")
+        expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+        assert output.dtype == expected_dtype
         assert np.all(output.shape == s.shape)
+
+        # Same test for the symmetric filter
         output = basis_function.filter_symmetric(s)
-        assert isinstance(output, np.ndarray)
-        assert np.issubdtype(output.dtype, np.floating)
+        assert array_api_compat.is_array_api_obj(output)
+        assert array_api_compat.array_namespace(output) is xp
+        # Integer input should be cast to the default floating point dtype of
+        # the array namespace while float32 and float64 input keeps its dtype.
+        assert xp.isdtype(output.dtype, kind="real floating")
+        expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+        assert output.dtype == expected_dtype
         assert np.all(output.shape == s.shape)
 
 
diff --git a/tests/test_spline_curves.py b/tests/test_spline_curves.py
index ead29df..8d22549 100644
--- a/tests/test_spline_curves.py
+++ b/tests/test_spline_curves.py
@@ -1,30 +1,54 @@
 import json
 import math
-import unittest.mock
 
+import array_api_compat
 import numpy as np
 import pytest
 import scipy
 import splinebox
 
 
-def test_check_control_points(non_hermite_spline_curve, coeff_gen):
+def test_get_control_points_xp_and_device_raises(non_hermite_spline_curve):
     spline = non_hermite_spline_curve
+    # Before setting the control points
     with pytest.raises(RuntimeError):
-        spline._check_control_points()
+        spline._get_control_points_xp_and_device()
+
+
+def test_get_control_points_xp_and_device(non_hermite_spline_curve, coeff_gen, xp, device):
+    spline = non_hermite_spline_curve
     spline.control_points = coeff_gen(spline.M, spline.basis_function.support, spline.closed)
-    spline._check_control_points()
+    returned_xp, returned_device = spline._get_control_points_xp_and_device()
+    assert returned_xp is xp
+    assert returned_device == device
 
 
-def test_check_control_points_and_tangents(hermite_spline_curve, coeff_gen):
-    spline = hermite_spline_curve
+def test_get_control_points_and_tangents_xp_and_device_raises(hermite_spline_curve, coeff_gen):
+    spline = hermite_spline_curve.copy()
+
+    # No control points and no tangents
     with pytest.raises(RuntimeError):
-        spline._check_control_points_and_tangents()
+        spline._get_control_points_and_tangents_xp_and_device()
+
+    # No tangents
     spline.control_points = coeff_gen(spline.M, spline.basis_function.support, spline.closed)
     with pytest.raises(RuntimeError):
-        spline._check_control_points_and_tangents()
+        spline._get_control_points_and_tangents_xp_and_device()
+
+    # No control points
+    spline = hermite_spline_curve.copy()
     spline.tangents = coeff_gen(spline.M, spline.basis_function.support, spline.closed)
-    spline._check_control_points_and_tangents()
+    with pytest.raises(RuntimeError):
+        spline._get_control_points_and_tangents_xp_and_device()
+
+
+def test_get_control_points_and_tangents_xp_and_device(hermite_spline_curve, coeff_gen, xp, device):
+    spline = hermite_spline_curve
+    spline.control_points = coeff_gen(spline.M, spline.basis_function.support, spline.closed)
+    spline.tangents = coeff_gen(spline.M, spline.basis_function.support, spline.closed)
+    returned_xp, returned_device = spline._get_control_points_and_tangents_xp_and_device()
+    assert returned_xp is xp
+    assert returned_device == device
 
 
 def test_minimum_number_of_knots(basis_function):
@@ -33,7 +57,7 @@ def test_minimum_number_of_knots(basis_function):
         splinebox.spline_curves.Spline(M, basis_function)
 
 
-def test_call(spline_curve, coeff_gen, derivative, call_positions, is_hermite_spline, is_interpolating):
+def test_call(spline_curve, coeff_gen, derivative, call_positions, is_hermite_spline, is_interpolating, xp):
     support = spline_curve.basis_function.support
     half_support = support / 2
     closed = spline_curve.closed
@@ -44,7 +68,8 @@ def test_call(spline_curve, coeff_gen, derivative, call_positions, is_hermite_sp
         spline_curve.tangents = coeff_gen(spline_curve.M, support, closed)
 
     if is_interpolating(spline_curve) and derivative == 0:
-        values = spline_curve(np.arange(spline_curve.M), derivative=derivative)
+        t = xp.arange(spline_curve.M, dtype=spline_curve.control_points.dtype)
+        values = spline_curve(t, derivative=derivative)
         expected = spline_curve.control_points
         if not closed:
             pad = math.ceil(half_support) - 1
@@ -55,11 +80,23 @@ def test_call(spline_curve, coeff_gen, derivative, call_positions, is_hermite_sp
     else:
         spline_curve(call_positions, derivative=derivative)
 
-    # Check that the presence of the coefficients and tangents was verified
-    if is_hermite_spline(spline_curve):
-        spline_curve._check_control_points_and_tangents.assert_called_once()
-    else:
-        spline_curve._check_control_points.assert_called_once()
+
+def test_raises_without_control_points(spline_curve, call_positions, translation_vector):
+    """
+    Methods that rely on the control points raise a RuntimeError
+    if they are called before the control points (and tangents) are set.
+    """
+    # The spline_curve fixture returns splines without control points set
+    with pytest.raises(RuntimeError):
+        spline_curve(call_positions)
+    with pytest.raises(RuntimeError):
+        spline_curve.arc_length(1.0)
+    with pytest.raises(RuntimeError):
+        spline_curve._control_points_centroid()
+    with pytest.raises(RuntimeError):
+        spline_curve.translate(translation_vector)
+    with pytest.raises(RuntimeError):
+        spline_curve.scale(1.5)
 
 
 def test_set_control_points(spline_curve, is_hermite_spline):
@@ -118,48 +155,48 @@ def test_hermite_spline_with_normal_basis(non_hermite_basis_function, M):
         splinebox.spline_curves.HermiteSpline(M, non_hermite_basis_function)
 
 
-def test_closed_splines(closed_spline_curve, derivative, coeff_gen, is_hermite_spline):
-    M = closed_spline_curve.M
-    support = closed_spline_curve.basis_function.support
-    closed_spline_curve.control_points = coeff_gen(M, support, closed=True)
-    if is_hermite_spline(closed_spline_curve):
-        closed_spline_curve.tangents = coeff_gen(M, support, closed=True)
-    assert np.allclose(
-        closed_spline_curve(0, derivative=derivative),
-        closed_spline_curve(M, derivative=derivative),
+def test_closed_splines(initialised_closed_spline_curve, derivative, xp):
+    spline = initialised_closed_spline_curve
+    assert xp.allclose(
+        spline(0, derivative=derivative),
+        spline(spline.M, derivative=derivative),
         equal_nan=True,
     )
 
 
-def test_draw(closed):
+def test_draw(closed, xp, device):
     spline = splinebox.spline_curves.Spline(M=4, basis_function=splinebox.basis_functions.B1(), closed=closed)
-    spline._check_control_points = unittest.mock.MagicMock()
-    spline.knots = np.array([[1, 1], [1, 2], [2, 2], [2, 1]])
+    spline.knots = xp.asarray([[1, 1], [1, 2], [2, 2], [2, 1]])
 
-    x = np.linspace(0, 3, 16)
-    y = np.linspace(0, 3, 16)
+    x = xp.linspace(0, 3, 16)
+    y = xp.linspace(0, 3, 16)
 
     # You can only draw closed splines
     if not closed:
         with pytest.raises(RuntimeError):
             spline.draw(x, y)
     else:
-        expected = np.zeros((len(x), len(y)))
-        expected[6:10, 6:10] = 1
-        expected[5, 5:11] = 0.5
-        expected[10, 5:11] = 0.5
-        expected[5:11, 5] = 0.5
-        expected[5:11, 10] = 0.5
+        expected = xp.zeros((len(x), len(y)))
+        if array_api_compat.is_jax_namespace(xp):
+            expected = expected.at[6:10, 6:10].set(1)
+            expected = expected.at[5, 5:11].set(0.5)
+            expected = expected.at[10, 5:11].set(0.5)
+            expected = expected.at[5:11, 5].set(0.5)
+            expected = expected.at[5:11, 10].set(0.5)
+        else:
+            expected[6:10, 6:10] = 1
+            expected[5, 5:11] = 0.5
+            expected[10, 5:11] = 0.5
+            expected[5:11, 5] = 0.5
+            expected[5:11, 10] = 0.5
 
         output = spline.draw(x, y)
-        assert np.allclose(output, expected)
 
-        # Check that the presence of the coefficients was verified
-        spline._check_control_points.assert_called()
+        assert xp.allclose(output, expected)
 
         # You can only draw splines in 2D
         spline = splinebox.spline_curves.Spline(M=4, basis_function=splinebox.basis_functions.B1(), closed=closed)
-        spline.knots = np.array([[1, 1, 0.5], [1, 2, 2], [2, 2, 1], [2, 1, 0]])
+        spline.knots = xp.asarray([[1, 1, 0.5], [1, 2, 2], [2, 2, 1], [2, 1, 0]])
         with pytest.raises(RuntimeError):
             spline.draw(x, y)
 
@@ -189,84 +226,70 @@ def test_is_inside():
         spline.is_inside(np.arange(10), np.arange(15))
 
 
-def test_arc_length():
+def test_arc_length(xp, device):
     # Create circular spline with radius sqrt(2)
     M = 10
     basis_function = splinebox.basis_functions.Exponential(M)
     spline = splinebox.spline_curves.Spline(M=M, basis_function=basis_function, closed=True)
-    spline._check_control_points = unittest.mock.MagicMock()
-    phi = np.linspace(0, 2 * np.pi, M + 1)[:-1]
-    knots = np.stack([np.cos(phi), np.sin(phi)], axis=-1)
+    phi = xp.linspace(0, 2 * np.pi, M + 1, device=device)[:-1]
+    knots = xp.stack([xp.cos(phi), xp.sin(phi)], axis=-1)
     spline.knots = knots
 
-    ts = np.linspace(0, M, 100)
-
-    # Reset mock in case it was called by any of the methods befor
-    spline._check_control_points.reset_mock()
+    ts = xp.linspace(0, M, 100, device=device)
 
-    arc_lengths = [spline.arc_length(t) for t in ts]
-    arc_lengths = np.array(arc_lengths)
+    arc_lengths = [spline.arc_length(t.item()) for t in ts]
+    arc_lengths = xp.asarray(arc_lengths, device=device)
 
-    # Check that the presence of the coefficients was verified
-    spline._check_control_points.assert_called()
-
-    expected = np.linspace(0, 2 * np.pi, 100)
-    assert np.allclose(arc_lengths, expected)
+    expected = xp.linspace(0, 2 * xp.pi, 100, device=device)
+    assert xp.allclose(arc_lengths, expected)
 
     # Check that it is working with a vector for stop
     permutation = np.random.permutation(len(ts))
+    permutation = xp.asarray(permutation, device=device)
     arc_lengths = spline.arc_length(ts[permutation])
-    assert np.allclose(arc_lengths, expected[permutation])
+    assert xp.allclose(arc_lengths, expected[permutation])
 
     # Check that is works with vectors for start and stop
-    arc_lengths = spline.arc_length(ts[permutation], np.zeros(len(ts)))
-    assert np.allclose(arc_lengths, expected[permutation])
+    arc_lengths = spline.arc_length(ts[permutation], xp.zeros(len(ts)))
+    assert xp.allclose(arc_lengths, expected[permutation])
 
 
-def test_arc_length_to_parameter():
+def test_arc_length_to_parameter(xp, device):
     # Create circular spline with radius sqrt(2)
     M = 100
     basis_function = splinebox.basis_functions.Exponential(M)
     spline = splinebox.spline_curves.Spline(M=M, basis_function=basis_function, closed=True)
-    spline._check_control_points = unittest.mock.MagicMock()
     knots = []
-    for t in np.linspace(0, 2 * np.pi, M + 1)[:-1]:
-        knots.append([np.cos(t), np.sin(t)])
-    knots = np.array(knots)
+    for t in xp.linspace(0, 2 * xp.pi, M + 1, device=device)[:-1]:
+        knots.append([xp.cos(t), xp.sin(t)])
+    knots = xp.asarray(knots, device=device)
     spline.knots = knots
 
-    ls = np.linspace(0, 2 * np.pi, 15)
-    expected = np.linspace(0, M, 15)  # ls / 2 / np.pi * M
+    ls = xp.linspace(0, 2 * xp.pi, 15, device=device)
+    expected = xp.linspace(0, M, 15, device=device)  # ls / 2 / xp.pi * M
     atol = 1e-5
 
-    # Reset mock in case it was called by any of the methods befor
-    spline._check_control_points.reset_mock()
-
     results = spline.arc_length_to_parameter(ls, atol=atol)
 
-    # Check that the presence of the coefficients was verified
-    spline._check_control_points.assert_called()
-
-    assert np.allclose(results, expected, atol=1e-3, rtol=0)
+    assert xp.allclose(results, expected, atol=1e-3, rtol=0)
 
     # Check that it works with an array of integer inputs
-    ls = np.arange(0, int(2 * np.pi), dtype=int)
-    expected = ls / 2 / np.pi * M
+    ls = xp.arange(0, int(2 * xp.pi), dtype=int, device=device)
+    expected = ls / 2 / xp.pi * M
     results = spline.arc_length_to_parameter(ls, atol=atol)
-    assert np.allclose(results, expected, atol=1e-3, rtol=0)
+    assert xp.allclose(results, expected, atol=1e-3, rtol=0)
     results = spline.arc_length_to_parameter(ls, atol=atol)
 
     # Create a sawtooth spline
     M = 7
     basis_function = splinebox.basis_functions.B1()
     spline = splinebox.spline_curves.Spline(M=M, basis_function=basis_function, closed=False)
-    spline._check_control_points = unittest.mock.MagicMock()
-    spline.knots = np.array([[0, 0], [1, 1], [1, 0], [2, 1], [2, 0], [3, 1], [3, 0]])
+    spline.knots = xp.asarray([[0, 0], [1, 1], [1, 0], [2, 1], [2, 0], [3, 1], [3, 0]], device=device)
 
-    rising_length = np.sqrt(2)
+    rising_length = math.sqrt(2)
     tooth_length = rising_length + 1
     total_length = 3 * tooth_length
-    ls = np.linspace(0, total_length, 15)
+    ls = xp.linspace(0, total_length, 15, device=device)
 
     def to_param(length):
         n_tooth = length // tooth_length
@@ -275,23 +298,17 @@ def to_param(length):
         # there are two knots per tooth
         return n_tooth * 2 + partial_tooth_param
 
-    expected = np.array(list(map(to_param, ls)))
+    expected = xp.asarray(list(map(to_param, ls)), device=device)
     atol = 1e-4
 
-    # Reset mock in case it was called by any of the methods befor
-    spline._check_control_points.reset_mock()
-
     results = spline.arc_length_to_parameter(ls, atol=atol)
 
-    # Check that the presence of the coefficients was verified
-    spline._check_control_points.assert_called()
-
-    assert np.allclose(results, expected, atol=atol)
+    assert xp.allclose(results, expected, atol=atol)
 
     # Check that is works for a single value
     expected = to_param(ls[10])
     result = spline.arc_length_to_parameter(ls[10], atol=atol)
-    assert np.isclose(result, expected, atol=atol)
+    assert xp.isclose(result, expected, atol=atol)
 
     # Check that length smaller than 0 and larger than the total length of the spline raise and error
     with pytest.raises(ValueError):
@@ -300,50 +317,60 @@ def to_param(length):
         spline.arc_length_to_parameter(total_length + 0.01)
 
 
-def test_translate(initialized_spline_curve, translation_vector):
-    spline = initialized_spline_curve
+def test_translate(initialised_spline_curve, translation_vector):
+    spline = initialised_spline_curve
     spline_copy = spline.copy()
-    spline_copy._check_control_points = unittest.mock.MagicMock()
     spline_copy.translate(translation_vector)
-    # Check that the presence of the coefficients was verified
-    spline_copy._check_control_points.assert_called()
     t = np.linspace(0, spline.M, 100) if spline.closed else np.linspace(0, spline.M - 1, 100)
     expected = spline(t) + translation_vector
 
     assert np.allclose(spline_copy(t), expected)
 
 
-def test_rotate(initialized_spline_curve, rotation_matrix, is_hermite_spline):
-    spline = initialized_spline_curve
-    spline_copy = spline.copy()
-    spline_copy._check_control_points = unittest.mock.MagicMock()
+def test_rotate(initialised_spline_curve, rotation_matrix, is_hermite_spline, xp, device):
+    spline = initialised_spline_curve
     if spline.ndim == 1:
         with pytest.raises(RuntimeError):
-            spline_copy.rotate(rotation_matrix)
+            spline.rotate(rotation_matrix)
     else:
+        spline_copy = spline.copy()
         spline_copy.rotate(rotation_matrix, centred=False)
-        # Check that the presence of the coefficients was verified
-        spline_copy._check_control_points.assert_called()
-        t = np.linspace(0, spline.M, 100) if spline.closed else np.linspace(0, spline.M - 1, 100)
+        t = (
+            xp.linspace(0, spline.M, 100, device=device)
+            if spline.closed
+            else xp.linspace(0, spline.M - 1, 100, device=device)
+        )
+        vals = spline(t)
+        expected = (rotation_matrix @ vals.T).T
+        assert xp.allclose(spline_copy(t), expected)
+
+        spline_copy = spline.copy()
+        spline_copy.rotate(rotation_matrix, centred=True)
+        t = (
+            xp.linspace(0, spline.M, 100, device=device)
+            if spline.closed
+            else xp.linspace(0, spline.M - 1, 100, device=device)
+        )
         vals = spline(t)
         expected = (rotation_matrix @ vals.T).T
-        assert np.allclose(spline_copy(t), expected)
+        assert xp.allclose(spline_copy(t), expected)
 
         spline_copy = spline.copy()
-        spline_copy._check_control_points = unittest.mock.MagicMock()
         spline_copy.rotate(rotation_matrix, centred=True)
-        # Check that the presence of the coefficients was verified
-        spline_copy._check_control_points.assert_called()
-        t = np.linspace(0, spline.M, 100) if spline.closed else np.linspace(0, spline.M - 1, 100)
+        t = (
+            xp.linspace(0, spline.M, 100, device=device)
+            if spline.closed
+            else xp.linspace(0, spline.M - 1, 100, device=device)
+        )
         vals = spline(t)
-        centring_vector = np.mean(spline.control_points, axis=0)
+        centring_vector = xp.mean(spline.control_points, axis=0)
         vals = vals - centring_vector
         expected = (rotation_matrix @ vals.T).T
         expected = expected + centring_vector
-        assert np.allclose(spline_copy(t), expected)
+        assert xp.allclose(spline_copy(t), expected)
 
 
-def test_fit(spline_curve, points, is_hermite_spline):
+def test_fit(spline_curve, points, is_hermite_spline, xp, device):
     # Bool indicating whether spline_curve is a hermite spline or not
     hermite = is_hermite_spline(spline_curve)
 
@@ -355,6 +382,14 @@ def test_fit(spline_curve, points, is_hermite_spline):
     else:
         spline_curve.fit(points)
 
+        assert array_api_compat.array_namespace(spline_curve.control_points) is xp
+        assert array_api_compat.device(spline_curve.control_points) == device
+
+        # To compare to the reference we use numpy
+        spline_curve.control_points = np.array(spline_curve.control_points)
+        if hermite:
+            spline_curve.tangents = np.array(spline_curve.tangents)
+
         # Keep a copy of the result of fit
         control_points_fit = spline_curve.control_points.copy()
         if hermite:
@@ -403,7 +438,7 @@ def difference_func(x, i):
             assert np.allclose(tangents_expected, tangents_fit)
 
 
-def test_fit_clamped(open_spline_curve, points, is_hermite_spline):
+def test_fit_clamped(open_spline_curve, points, is_hermite_spline, xp):
     spline = open_spline_curve
 
     # Bool indicating whether spline is a hermite spline or not
@@ -423,8 +458,8 @@ def test_fit_clamped(open_spline_curve, points, is_hermite_spline):
         spline.fit(points, "clamped")
 
         # Check that the boundary conditions are fulfilled
-        assert np.allclose(spline(0, derivative=1), 0)
-        assert np.allclose(spline(spline.M - 1, derivative=1), 0)
+        assert xp.allclose(spline(0, derivative=1), xp.asarray([0], dtype=points.dtype))
+        assert xp.allclose(spline(spline.M - 1, derivative=1), xp.asarray([0], dtype=points.dtype))
 
 
 def test_fit_natural(open_spline_curve, points, is_hermite_spline):
@@ -456,7 +491,7 @@ def test_fit_natural(open_spline_curve, points, is_hermite_spline):
         assert np.allclose(spline(spline.M - 1, derivative=2), 0)
 
 
-def test_knots(spline_curve, knot_gen, is_hermite_spline, request):
+def test_knots(spline_curve, knot_gen, is_hermite_spline, request, xp, dtype):
     if isinstance(spline_curve.basis_function, splinebox.basis_functions.B2):
         # The filter_symmetric and filter_periodic are not implemented for B2
         request.node.add_marker(pytest.mark.xfail)
@@ -471,26 +506,34 @@ def test_knots(spline_curve, knot_gen, is_hermite_spline, request):
         # Use finite differences as tangents
         if spline_curve.closed:
             if knots.ndim == 1:
-                diff = np.diff(knots, axis=0, append=knots[0])
+                diff = xp.diff(knots, axis=0, append=knots[0])
             else:
-                diff = np.diff(knots, axis=0, append=knots[0][np.newaxis, :])
+                diff = xp.diff(knots, axis=0, append=knots[0][np.newaxis, :])
             tangents = (diff + np.roll(diff, 1)) / 2
         else:
-            diff = np.diff(knots, axis=0)
+            diff = xp.diff(knots, axis=0)
             pad = math.ceil(spline_curve.basis_function.support / 2) - 1
             if diff.ndim == 1:
                 diff = diff[:, np.newaxis]
-            diff = np.pad(diff, ((pad + 1, pad + 1), (0, 0)), mode="constant", constant_values=(0,))
+            padding = xp.tile(xp.zeros(diff.shape[1], dtype=diff.dtype), (pad + 1, 1))
+            diff = xp.concat((padding, diff, padding), axis=0)
             tangents = (diff[:-1] + diff[1:]) / 2
-        spline_curve.tangents = tangents
+        spline_curve.tangents = xp.astype(tangents, dtype)
+
+    if dtype is int:
+        # For comparison the knots have to be of the same type as the returned values.
+        # For int the default floating point type is returned.
+        knots = xp.astype(knots, xp.asarray(1.0).dtype)
+
     # This assert makes sense because the knots are not saved in
     # the spline class but are converted to coefficients which are saved.
+    atol = 1e-5 if knots.dtype == xp.float32 else 1e-8
     if spline_curve.padding_function is None and not spline_curve.closed and pad != 0:
         # We only compare the actual knots because padded knots will not
         # be the same because of the filtering to find the control points.
-        assert np.allclose(knots[pad:-pad], spline_curve.knots[pad:-pad])
+        assert xp.allclose(knots[pad:-pad], spline_curve.knots[pad:-pad], atol=atol)
     else:
-        assert np.allclose(knots, spline_curve.knots)
+        assert xp.allclose(knots, spline_curve.knots, atol=atol)
 
 
 def test_number_of_knots(spline_curve, knot_gen):
@@ -513,7 +556,7 @@ def test_number_of_knots(spline_curve, knot_gen):
             spline_curve.knots = knots
 
 
-def test_centroid(spline_curve, coeff_gen, is_hermite_spline):
+def test_centroid(spline_curve, coeff_gen, is_hermite_spline, xp, dtype):
     support = spline_curve.basis_function.support
     closed = spline_curve.closed
     M = spline_curve.M
@@ -524,14 +567,13 @@ def test_centroid(spline_curve, coeff_gen, is_hermite_spline):
         tangents = coeff_gen(spline_curve.M, support, closed)
         spline_curve.tangents = tangents
 
-    expected = np.mean(control_points, axis=0)
-    assert np.allclose(spline_curve._control_points_centroid(), expected)
+    output_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    expected = xp.mean(control_points, axis=0, dtype=output_dtype)
+    assert xp.allclose(spline_curve._control_points_centroid(), expected)
 
-    spline_curve._check_control_points.assert_called()
 
-
-def test_scale(initialized_spline_curve, is_hermite_spline):
-    spline = initialized_spline_curve
+def test_scale(initialised_spline_curve, is_hermite_spline, xp):
+    spline = initialised_spline_curve
     spline_copy = spline.copy()
 
     factor = 1.5
@@ -544,30 +586,18 @@ def test_scale(initialized_spline_curve, is_hermite_spline):
     if is_hermite_spline(spline_copy):
         spline_copy.tangents = spline_copy.tangents * factor
 
-    t = np.linspace(0, spline.M, 100) if spline.closed else np.linspace(0, spline.M - 1, 100)
-
-    # Reset mocks in case they were called by another method already
-    if is_hermite_spline(spline_copy):
-        spline._check_control_points_and_tangents.reset_mock()
-    else:
-        spline._check_control_points.reset_mock()
+    t = xp.linspace(0, spline.M, 100) if spline.closed else xp.linspace(0, spline.M - 1, 100)
 
     spline.scale(factor)
 
-    # Check that coefficients and tangents were checked
-    if is_hermite_spline(spline_copy):
-        spline._check_control_points_and_tangents.assert_called()
-    else:
-        spline._check_control_points.assert_called()
+    assert xp.allclose(spline(t), spline_copy(t))
 
-    assert np.allclose(spline(t), spline_copy(t))
 
-
-def test_curvilinear_reparametrisation_energy():
+def test_curvilinear_reparametrisation_energy(as_xp_array, xp):
     M = 4
     spline = splinebox.spline_curves.Spline(M, splinebox.B1())
 
-    spline.knots = np.array([[1, 0], [2, 0], [3, 0], [4, 0]])
+    spline.knots = as_xp_array([[1, 0], [2, 0], [3, 0], [4, 0]])
 
     # Because of the spacing of the knots the derivative should be constant so we can just compute one value in the middle of the spline.
     derivative_val = spline((spline.M - 1) / 2, derivative=1)
@@ -576,145 +606,156 @@ def test_curvilinear_reparametrisation_energy():
     # We can replace the integration in the definition of the curvilinear reparametrisation energy with a multiplication because the derivative is constant.
     length = spline.arc_length()
     c = (length / spline.M) ** 2
-    e_curv = (np.linalg.norm(derivative_val) ** 2 - c) ** 2
+    e_curv = (xp.linalg.vector_norm(derivative_val) ** 2 - c) ** 2
     expected_val = e_curv * (spline.M - 1) / length**4
 
     val = spline.curvilinear_reparametrisation_energy()
 
-    assert np.isclose(val, expected_val)
+    assert xp.isclose(val, expected_val)
 
 
-def test_curvilinear_reparametrisation_energy_user_c():
+def test_curvilinear_reparametrisation_energy_user_c(as_xp_array, xp):
     M = 4
     spline = splinebox.spline_curves.Spline(M, splinebox.B1())
 
-    spline.knots = np.array([[1, 0], [2, 0], [3, 0], [4, 0]])
+    spline.knots = as_xp_array([[1, 0], [2, 0], [3, 0], [4, 0]])
 
     # Because of the spacing of the knots the derivative should be constant so we can just compute one value in the middle of the spline.
     derivative_val = spline((spline.M - 1) / 2, derivative=1)
 
     # For a user provided c the unnormalised integral (equation 25) is returned.
     c = 1.5
-    expected_val = (np.linalg.norm(derivative_val) ** 2 - c) ** 2 * (spline.M - 1)
+    expected_val = (xp.linalg.vector_norm(derivative_val) ** 2 - c) ** 2 * (spline.M - 1)
 
     val = spline.curvilinear_reparametrisation_energy(c=c)
 
-    assert np.isclose(val, expected_val)
+    assert xp.isclose(val, expected_val)
 
 
-def test_curvilinear_reparametrisation_energy_translation(initialized_spline_curve, translation_vector):
+def test_curvilinear_reparametrisation_energy_translation(initialised_spline_curve, translation_vector, xp):
     """
     Test if the curvilinear reparametrisation energy is invariant to translation.
     """
     atol = 1e-3
-    spline = initialized_spline_curve
+    spline = initialised_spline_curve
     expected = spline.curvilinear_reparametrisation_energy(atol=atol)
     spline.translate(translation_vector)
     val = spline.curvilinear_reparametrisation_energy(atol=atol)
-    assert np.isclose(val, expected, atol=atol)
+    assert xp.isclose(val, expected, atol=atol)
 
 
-def test_curvilinear_reparametrisation_energy_scale_invariance(initialized_spline_curve):
+def test_curvilinear_reparametrisation_energy_scale_invariance(initialised_spline_curve, xp):
     """
     The curvilinear reparametrisation energy should not change when the entire spline is scaled.
     """
-    spline = initialized_spline_curve
+    spline = initialised_spline_curve
     expected = spline.curvilinear_reparametrisation_energy()
     spline.scale(0.1)
     val = spline.curvilinear_reparametrisation_energy()
-    assert np.isclose(val, expected)
+    assert xp.isclose(val, expected)
 
 
-def test_curvature():
+def test_curvature(as_xp_array, xp):
     # Create a circular spline
     M = 4
     spline = splinebox.spline_curves.Spline(M=M, basis_function=splinebox.basis_functions.Exponential(M), closed=True)
-    t = np.linspace(0, M, 100)
+    t = xp.linspace(0, M, 100)
     for r in (1, 1.5, 2):
         # 2D
-        spline.knots = np.array([[0, r], [r, 0], [0, -r], [-r, 0]])
+        spline.knots = as_xp_array([[0, r], [r, 0], [0, -r], [-r, 0]])
         curvature = spline.curvature(t)
-        assert np.allclose(curvature, 1 / r)
+        assert xp.allclose(curvature, 1 / r)
 
         # If we reverse the direction, the curvature should change sign
-        spline.knots = np.array([[0, r], [-r, 0], [0, -r], [r, 0]])
+        spline.knots = as_xp_array([[0, r], [-r, 0], [0, -r], [r, 0]])
         curvature = spline.curvature(t)
-        assert np.allclose(curvature, -1 / r)
+        assert xp.allclose(curvature, -1 / r)
 
         # 3D
-        spline.knots = np.array([[0, 0, r], [0, r, 0], [0, 0, -r], [0, -r, 0]])
+        spline.knots = as_xp_array([[0, 0, r], [0, r, 0], [0, 0, -r], [0, -r, 0]])
         curvature = spline.curvature(t)
-        assert np.allclose(curvature, 1 / r)
+        assert xp.allclose(curvature, 1 / r)
 
     # 1D
     spline = splinebox.spline_curves.Spline(M=M, basis_function=splinebox.basis_functions.B3(), closed=False)
-    spline.knots = np.array([[1], [5], [4], [2]])
+    spline.knots = as_xp_array([[1], [5], [4], [2]])
     expected = spline(t, derivative=2) / (1 + spline(t, derivative=1) ** 2) ** (3 / 2)
-    expected = np.squeeze(expected)
+    expected = xp.squeeze(expected, axis=1)
     result = spline.curvature(t)
-    assert np.allclose(result, expected)
+    assert xp.allclose(result, expected)
 
 
-def test_normal_2D():
+def test_normal_2D(as_xp_array, xp):
     # Create a circular spline
     M = 4
     spline = splinebox.spline_curves.Spline(M=M, basis_function=splinebox.basis_functions.Exponential(M), closed=True)
-    t = np.linspace(0, M, 100)
+    t = xp.linspace(0, M, 100)
 
-    theta = np.linspace(0, 2 * np.pi, 100)
-    expected = np.stack([np.sin(theta), np.cos(theta)], axis=-1)
+    theta = xp.linspace(0, 2 * xp.pi, 100)
+    expected = xp.stack([xp.sin(theta), xp.cos(theta)], axis=-1)
     for r in (1, 1.5, 2):
-        spline.knots = np.array([[0, r], [r, 0], [0, -r], [-r, 0]])
+        spline.knots = as_xp_array([[0, r], [r, 0], [0, -r], [-r, 0]])
         normals = spline.normal(t)
-        assert np.allclose(normals, expected)
+        assert xp.allclose(normals, expected)
 
         # Reversing the direction change the direction of the normal vector
-        spline.knots = np.array([[0, r], [-r, 0], [0, -r], [r, 0]])
+        spline.knots = as_xp_array([[0, r], [-r, 0], [0, -r], [r, 0]])
         normals = spline.normal(t[::-1])
-        assert np.allclose(normals, -1 * expected)
+        assert xp.allclose(normals, -1 * expected)
 
         # Check that they are normal vectors
-        assert np.allclose(np.linalg.norm(normals, axis=1), np.ones(len(t)))
+        assert xp.allclose(xp.linalg.vector_norm(normals, axis=1), xp.ones(len(t)))
 
 
-def test_normal_3D():
+def test_normal_3D(as_xp_array, xp):
     M = 4
     spline = splinebox.spline_curves.Spline(M=M, basis_function=splinebox.basis_functions.B3(), closed=False)
-    t = np.linspace(0, M - 1, 100)
+    t = xp.linspace(0, M - 1, 100)
 
     # This is a straight line
-    spline.control_points = np.stack([np.arange(M + 2), np.arange(M + 2), np.ones(M + 2)]).T
-    normals = spline.normal(t, frame="bishop", initial_vector=np.array([0.0, 0.0, 1.0]))
+    spline.control_points = xp.stack([xp.arange(M + 2), xp.arange(M + 2), xp.ones(M + 2)]).T
+    normals = spline.normal(t, frame="bishop", initial_vector=as_xp_array([0.0, 0.0, 1.0]))
     # check that the normals don't change
-    assert np.isclose(np.sum(np.diff(normals, axis=0)), 0)
+    assert xp.isclose(xp.sum(xp.diff(normals, axis=0)), 0)
 
-    spline.knots = np.array([[1, 0, 0], [0, 1, 1], [-1, 0, 2], [0, -1, 3]])
+    spline.knots = as_xp_array([[1, 0, 0], [0, 1, 1], [-1, 0, 2], [0, -1, 3]])
     t = t[1:-1]
     normals = spline.normal(t, frame="frenet")
     tangents = spline(t, derivative=1)
     # Check that both normals are orthogonal t the tangents (dot product)
-    dot = np.sum(normals * tangents[:, np.newaxis], axis=-1)
-    assert np.allclose(dot, 0)
+    dot = xp.sum(normals * tangents[:, xp.newaxis], axis=-1)
+    assert xp.allclose(dot, 0)
     # Check that the normals are orthogonal to each other
-    dot = np.sum(normals[:, 0] * normals[:, 1])
-    assert np.isclose(dot, 0)
+    dot = xp.sum(normals[:, 0] * normals[:, 1])
+    assert xp.isclose(dot, 0)
 
 
-def test_repr(initialized_spline_curve, is_hermite_spline):
-    spline = initialized_spline_curve
-    # The rounding is necessary because the precision is limited in the string representation
-    spline.control_points = np.round(spline.control_points, decimals=np.get_printoptions()["precision"])
-    if is_hermite_spline(spline):
-        spline.tangents = np.round(spline.tangents, decimals=np.get_printoptions()["precision"])
+def test_repr(initialised_spline_curve, is_hermite_spline, xp, dtype):
+    spline = initialised_spline_curve
+
+    if dtype is not int:
+        # The rounding is necessary because the precision is limited in the string representation
+        xp.set_printoptions(precision=16)
+        spline.control_points = xp.round(spline.control_points, decimals=1)
+        if is_hermite_spline(spline):
+            spline.tangents = xp.round(spline.tangents, decimals=1)
+    if array_api_compat.is_jax_namespace(xp):
+        pass
     assert eval(repr(spline)) == spline
 
 
 def test_eq_spline(spline_curve, basis_function, closed, coeff_gen, is_hermite_spline):
     spline = spline_curve
     basis_function = basis_function
+    spline.control_points = coeff_gen(spline.M, spline.basis_function.support, spline.closed)
     if is_hermite_spline(spline):
+        spline.tangents = coeff_gen(spline.M, spline.basis_function.support, spline.closed)
         other = splinebox.spline_curves.HermiteSpline(
-            M=spline.M, basis_function=basis_function, closed=closed, control_points=spline.control_points
+            M=spline.M,
+            basis_function=basis_function,
+            closed=closed,
+            control_points=spline.control_points,
+            tangents=spline.tangents,
         )
     else:
         other = splinebox.spline_curves.Spline(
@@ -726,8 +767,6 @@ def test_eq_spline(spline_curve, basis_function, closed, coeff_gen, is_hermite_s
         assert spline == other
 
         spline.control_points = coeff_gen(spline.M, spline.basis_function.support, spline.closed)
-        if is_hermite_spline(spline):
-            spline.tangents = coeff_gen(spline.M, spline.basis_function.support, spline.closed)
         assert spline != other
 
         other.control_points = spline.control_points
@@ -748,15 +787,15 @@ def test_eq_spline(spline_curve, basis_function, closed, coeff_gen, is_hermite_s
     assert spline != other
 
 
-def test_str(initialized_spline_curve):
-    assert isinstance(str(initialized_spline_curve), str)
+def test_str(initialised_spline_curve):
+    assert isinstance(str(initialised_spline_curve), str)
 
 
-def test_to_json_and_from_json(initialized_spline_curve, tmpdir, is_hermite_spline):
-    spline = initialized_spline_curve
+def test_to_json_and_from_json(initialised_spline_curve, tmpdir, is_hermite_spline):
+    spline = initialised_spline_curve
     path = tmpdir / "spline.json"
     spline.to_json(path)
-    if is_hermite_spline(initialized_spline_curve):
+    if is_hermite_spline(initialised_spline_curve):
         loaded_spline_curve = splinebox.spline_curves.HermiteSpline.from_json(path)
     else:
         loaded_spline_curve = splinebox.spline_curves.Spline.from_json(path)
@@ -824,47 +863,65 @@ def test_saving_and_loading_of_multiple_splines(is_hermite_basis_function, tmpdi
         assert spline == loaded_spline
 
 
-def test_distance(initialized_spline_curve):
-    spline = initialized_spline_curve
+def test_distance(initialised_spline_curve, xp, device, request):
+    spline = initialised_spline_curve
+
+    if not array_api_compat.is_numpy_namespace(xp) and spline.ndim != 1:
+        # Currently only implemented for numpy
+        request.node.add_marker(pytest.mark.xfail)
 
     if spline.ndim == 1:
         with pytest.raises(RuntimeError):
-            spline.distance(np.array((1.0,)))
+            spline.distance(xp.asarray((1.0,), device=device))
         return
 
-    point = np.array([1] * spline.control_points.shape[1])
+    point = xp.asarray([1] * spline.control_points.shape[1], device=device)
     distance = spline.distance(point)
 
     # Check that a set of equaly spaced points on the spline are at least as far away as the distance
-    t = np.linspace(0, spline.M if spline.closed else spline.M - 1, spline.M * 10)
+    t = xp.linspace(0, spline.M if spline.closed else spline.M - 1, spline.M * 10, device=device)
     points_on_spline = spline(t)
-    distances = np.linalg.norm(points_on_spline - point[np.newaxis], axis=-1)
-    assert np.all(distances >= distance)
+    distances = xp.linalg.vector_norm(points_on_spline - point[xp.newaxis], axis=-1)
+    assert xp.all(distances >= distance)
 
     # Check that the returned t corresponds to the correct distance
     distance, t = spline.distance(point, return_t=True)
     point_on_spline = spline(t)
-    assert np.isclose(distance, np.linalg.norm(point - point_on_spline))
+    assert xp.isclose(distance, xp.linalg.vector_norm(point - point_on_spline))
 
 
-def test_distance_multiple_points(initialized_spline_curve):
-    spline = initialized_spline_curve
+def test_distance_multiple_points(initialised_spline_curve, as_xp_array, xp, request):
+    spline = initialised_spline_curve
+
+    if not array_api_compat.is_numpy_namespace(xp) and spline.ndim != 1:
+        # Currently only implemented for numpy
+        request.node.add_marker(pytest.mark.xfail)
 
     if spline.ndim == 1:
         return
 
     np.random.seed(45)
     points = np.random.rand(10, spline.control_points.shape[1]) * 10
+    points = as_xp_array(points)
 
     distances, t = spline.distance(points, return_t=True)
 
     points_on_spline = spline(t)
-    assert np.allclose(distances, np.linalg.norm(points - points_on_spline, axis=1))
+    assert xp.allclose(distances, xp.linalg.vector_norm(points - points_on_spline, axis=1))
+
+
+def test_moving_frame(initialised_spline_curve, xp, device, dtype):
+    spline = initialised_spline_curve
+    t = xp.linspace(0, spline.M if spline.closed else spline.M - 1, spline.M * 10, dtype=dtype, device=device)
+
+    # The dtype actually used. If dtype is int,
+    # the default float precision, otherwise the just dtype.
+    effective_dtype = xp.asarray(1.0, device=device).dtype if dtype is int else dtype
 
+    atol = 1e-5 if effective_dtype == xp.float32 else 1e-8
 
-def test_moving_frame(initialized_spline_curve):
-    spline = initialized_spline_curve
-    t = np.linspace(0, spline.M if spline.closed else spline.M - 1, spline.M * 10)
+    # pytorch requires a tensor for comparison
+    zero = xp.asarray(0.0, dtype=effective_dtype, device=device)
 
     if spline.ndim != 3:
         with pytest.raises(RuntimeError):
@@ -878,25 +935,30 @@ def test_moving_frame(initialized_spline_curve):
             with pytest.raises(RuntimeError):
                 frame = spline.moving_frame(t, method="frenet")
         else:
-            # The second derivative of B1 splines at the knots is not defined so we
+            # The second derivative of CatmullRom splines at the knots is not defined so we
             # remove all of the integers from t
             t_frenet = t[(t % 1) != 0] if isinstance(spline.basis_function, splinebox.basis_functions.CatmullRom) else t
+            # When dtype is int no t values are left
+            if len(t_frenet) == 0:
+                return
 
             frame = spline.moving_frame(t_frenet, method="frenet")
 
             # Check that you can pass a float value
-            assert np.allclose(frame[5], spline.moving_frame(t_frenet[5], method="frenet"))
+            assert xp.allclose(frame[5], spline.moving_frame(t_frenet[5], method="frenet"), atol=atol)
 
             # Check that they are unit vectors
-            assert np.allclose(np.linalg.norm(frame, axis=-1), 1)
+            assert xp.allclose(
+                xp.linalg.vector_norm(frame, axis=-1), xp.asarray(1.0, dtype=frame.dtype, device=device), atol=atol
+            )
 
             # Check if the frame forms an orthogonal basis using the dot product
-            assert np.allclose(np.sum(frame[:, 0] * frame[:, 1], axis=-1), 0)
-            assert np.allclose(np.sum(frame[:, 0] * frame[:, 2], axis=-1), 0)
-            assert np.allclose(np.sum(frame[:, 1] * frame[:, 2], axis=-1), 0)
+            assert xp.allclose(xp.sum(frame[:, 0] * frame[:, 1], axis=-1), zero, atol=atol)
+            assert xp.allclose(xp.sum(frame[:, 0] * frame[:, 2], axis=-1), zero, atol=atol)
+            assert xp.allclose(xp.sum(frame[:, 1] * frame[:, 2], axis=-1), zero, atol=atol)
 
             # Check that T is parallel to the first derivative
-            assert np.allclose(np.cross(frame[:, 0], spline(t_frenet, derivative=1)), 0)
+            assert xp.allclose(xp.cross(frame[:, 0], spline(t_frenet, derivative=1)), zero, atol=atol)
 
             # N = dT/ds and T = dr/ds so N should be parallel to d^2r/ds^2 - proj_T(d^2/ds^2)
             # Removing the projection on T is equivalent to removing the part of the
@@ -905,8 +967,10 @@ def test_moving_frame(initialized_spline_curve):
             # unit vector T, which has constant length.
             deriv2 = spline(t_frenet, derivative=2)
             # use Gram Schmidt to remove the scaling component
-            deriv2 -= np.sum(deriv2 * frame[:, 0], axis=-1)[:, np.newaxis] * frame[:, 0]
-            assert np.allclose(np.cross(frame[:, 1], deriv2), 0)
+            deriv2 -= xp.sum(deriv2 * frame[:, 0], axis=-1)[:, xp.newaxis] * frame[:, 0]
+            assert xp.allclose(
+                xp.cross(frame[:, 1], deriv2), zero, atol=atol * float(xp.max(xp.linalg.vector_norm(deriv2, axis=-1)))
+            )
 
         ########
         # Bishop
@@ -915,38 +979,49 @@ def test_moving_frame(initialized_spline_curve):
             # The derivative of B1 splines at the knots is not defined so we
             # remove all of the integers from t
             t = t[(t % 1) != 0]
+            # When dtype is int no t values are left
+            if len(t) == 0:
+                return
 
         with pytest.raises(ValueError):
             # Non-orthogonal initial vector
-            frame = spline.moving_frame(t, method="bishop", initial_vector=np.ones(3))
+            frame = spline.moving_frame(t, method="bishop", initial_vector=xp.ones(3, dtype=dtype, device=device))
 
-        initial_vector = np.zeros(3)
-        tangent = spline(t[0], derivative=1)
-        initial_vector[1] = -tangent[2]
-        initial_vector[2] = tangent[1]
+        tangent = spline(t[0].item(), derivative=1)
+        initial_vector = xp.zeros(3, dtype=tangent.dtype, device=device)
+        if array_api_compat.is_jax_namespace(xp):
+            initial_vector = initial_vector.at[1].set(-tangent[2])
+            initial_vector = initial_vector.at[2].set(tangent[1])
+        else:
+            initial_vector[1] = -tangent[2]
+            initial_vector[2] = tangent[1]
         frame = spline.moving_frame(t, method="bishop", initial_vector=initial_vector)
 
         # Check that you can pass a float value
-        assert np.allclose(frame[5], spline.moving_frame(t[5], method="bishop", initial_vector=frame[5, 1]))
+        assert xp.allclose(frame[5], spline.moving_frame(t[5], method="bishop", initial_vector=frame[5, 1]), atol=atol)
 
         # Check that they are unit vectors
-        assert np.allclose(np.linalg.norm(frame, axis=-1), 1)
+        assert xp.allclose(
+            xp.linalg.vector_norm(frame, axis=-1), xp.asarray(1.0, dtype=frame.dtype, device=device), atol=atol
+        )
 
         # Check if the frame forms an orthogonal basis using the dot product
-        assert np.allclose(np.sum(frame[:, 0] * frame[:, 1], axis=-1), 0)
-        assert np.allclose(np.sum(frame[:, 0] * frame[:, 2], axis=-1), 0)
-        assert np.allclose(np.sum(frame[:, 1] * frame[:, 2], axis=-1), 0)
+        # Since the Bishop frame is build iteratively, it accumulates more error.
+        atol_bishop = 1e-3 if effective_dtype == xp.float32 else 1e-8
+        assert xp.allclose(xp.sum(frame[:, 0] * frame[:, 1], axis=-1), zero, atol=atol_bishop)
+        assert xp.allclose(xp.sum(frame[:, 0] * frame[:, 2], axis=-1), zero, atol=atol_bishop)
+        assert xp.allclose(xp.sum(frame[:, 1] * frame[:, 2], axis=-1), zero, atol=atol_bishop)
 
         # Check that T is parallel to the first derivative
-        assert np.allclose(np.cross(frame[:, 0], spline(t, derivative=1)), 0)
+        assert xp.allclose(xp.cross(frame[:, 0], spline(t, derivative=1)), zero, atol=atol)
 
 
-def test_moving_frame_raises():
+def test_moving_frame_raises(as_xp_array, xp):
     M = 6
     spline = splinebox.spline_curves.Spline(M=M, basis_function=splinebox.basis_functions.B3(), closed=False)
-    t = np.linspace(0, M - 1, 100)
+    t = xp.linspace(0, M - 1, 100)
 
-    spline.knots = np.random.rand(6, 3)
+    spline.knots = as_xp_array(np.random.rand(6, 3))
 
     # Wrong method
     with pytest.raises(ValueError):
@@ -957,32 +1032,36 @@ def test_moving_frame_raises():
         spline.moving_frame(t, method="frenet")
 
     # Straight segment
-    spline.control_points = np.stack([np.arange(M + 2), np.zeros(M + 2), np.zeros(M + 2)], axis=-1)
+    spline.control_points = xp.stack([xp.arange(M + 2), xp.zeros(M + 2), xp.zeros(M + 2)], axis=-1)
     with pytest.raises(RuntimeError):
         spline.moving_frame(t, method="frenet")
 
 
-def test_moving_frame_unsorted(frame):
+def test_moving_frame_unsorted(frame, as_xp_array, xp):
     np.random.seed(499)
 
     spline = splinebox.spline_curves.Spline(M=7, basis_function=splinebox.basis_functions.B3(), closed=False)
-    spline.knots = np.random.rand(7, 3)
+    spline.knots = as_xp_array(np.random.rand(7, 3))
 
-    t = np.random.rand(100) * (spline.M - 1)
+    t = as_xp_array(np.random.rand(100)) * (spline.M - 1)
 
-    sort_indices = np.argsort(t)
+    sort_indices = xp.argsort(t)
     sorted_t = t[sort_indices]
 
-    initial_vector = np.cross(spline(t[0], derivative=1), spline(sorted_t[0], derivative=1))
+    initial_vector = xp.cross(spline(t[0], derivative=1), spline(sorted_t[0], derivative=1))
 
     sorted_frame = spline.moving_frame(sorted_t, method=frame, initial_vector=initial_vector)
     unsorted_frame = spline.moving_frame(t, method=frame, initial_vector=initial_vector)
 
-    assert np.allclose(sorted_frame[np.argsort(sort_indices)], unsorted_frame)
+    assert xp.allclose(sorted_frame[xp.argsort(sort_indices)], unsorted_frame)
+
 
+def test_mesh(initialised_spline_curve, radius, step_t, step_angle, mesh_type, cap_ends, frame, xp, request):
+    spline = initialised_spline_curve
 
-def test_mesh(initialized_spline_curve, radius, step_t, step_angle, mesh_type, cap_ends, frame):
-    spline = initialized_spline_curve
+    if not array_api_compat.is_numpy_namespace(xp) and spline.ndim == 3:
+        # Currently only implemented for numpy
+        request.node.add_marker(pytest.mark.xfail)
 
     if spline.ndim != 3:
         # mesh is only implemented for splines in 3D
@@ -1002,7 +1081,7 @@ def test_mesh(initialized_spline_curve, radius, step_t, step_angle, mesh_type, c
         )
 
         # Check that all points are part of at least one element
-        assert np.all(np.arange(len(points)) == np.sort(np.unique(connectivity)))
+        assert xp.all(xp.arange(len(points)) == xp.sort(xp.unique(connectivity)))
 
         if mesh_type == "surface" and (cap_ends is not None or spline.closed):
             # Compute the Euler characteristic to check that the mesh is closed
@@ -1031,29 +1110,29 @@ def test_mesh(initialized_spline_curve, radius, step_t, step_angle, mesh_type, c
                 assert V + F - E == 2
 
 
-def test_mesh_cap_ends_false_alias_matches_none():
+def test_mesh_cap_ends_false_alias_matches_none(as_xp_array, xp):
     spline = splinebox.spline_curves.Spline(M=4, basis_function=splinebox.basis_functions.B3(), closed=False)
-    spline.knots = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.5, 0.0], [3.0, 0.5, 0.5]])
+    spline.knots = as_xp_array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.5, 0.0], [3.0, 0.5, 0.5]])
 
     points_none, connectivity_none = spline.mesh(radius=0.2, step_t=0.5, step_angle=30, cap_ends=None)
     with pytest.warns(DeprecationWarning, match="cap_ends=False"):
         points_false, connectivity_false = spline.mesh(radius=0.2, step_t=0.5, step_angle=30, cap_ends=False)
 
-    assert np.allclose(points_none, points_false)
-    assert np.array_equal(connectivity_none, connectivity_false)
+    assert xp.allclose(points_none, points_false)
+    assert xp.array_equal(connectivity_none, connectivity_false)
 
 
-def test_mesh_cap_ends_true_raises_value_error():
+def test_mesh_cap_ends_true_raises_value_error(as_xp_array):
     spline = splinebox.spline_curves.Spline(M=4, basis_function=splinebox.basis_functions.B3(), closed=False)
-    spline.knots = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.5, 0.0], [3.0, 0.5, 0.5]])
+    spline.knots = as_xp_array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.5, 0.0], [3.0, 0.5, 0.5]])
 
     with pytest.raises(ValueError, match="cap_ends='flat'"):
         spline.mesh(radius=0.2, step_t=0.5, step_angle=30, cap_ends=True)
 
 
-def test_spherical_caps_add_more_geometry_than_flat_caps():
+def test_spherical_caps_add_more_geometry_than_flat_caps(as_xp_array):
     spline = splinebox.spline_curves.Spline(M=4, basis_function=splinebox.basis_functions.B3(), closed=False)
-    spline.knots = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.5, 0.0], [3.0, 0.5, 0.5]])
+    spline.knots = as_xp_array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.5, 0.0], [3.0, 0.5, 0.5]])
 
     flat_points, flat_connectivity = spline.mesh(radius=0.2, step_t=0.5, step_angle=30, cap_ends="flat")
     sphere_points, sphere_connectivity = spline.mesh(radius=0.2, step_t=0.5, step_angle=30, cap_ends="sphere")
@@ -1062,14 +1141,14 @@ def test_spherical_caps_add_more_geometry_than_flat_caps():
     assert len(sphere_connectivity) > len(flat_connectivity)
 
 
-def test_flat_caps_use_exact_open_spline_endpoints_when_step_t_does_not_divide_range():
+def test_flat_caps_use_exact_open_spline_endpoints_when_step_t_does_not_divide_range(as_xp_array, xp):
     spline = splinebox.spline_curves.Spline(M=4, basis_function=splinebox.basis_functions.B3(), closed=False)
-    spline.knots = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.5, 0.0], [3.0, 0.5, 0.5]])
+    spline.knots = as_xp_array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.5, 0.0], [3.0, 0.5, 0.5]])
 
     points, _ = spline.mesh(radius=0.2, step_t=0.4, step_angle=30, cap_ends="flat")
 
-    assert np.allclose(points[-2], spline(0))
-    assert np.allclose(points[-1], spline(spline.M - 1))
+    assert xp.allclose(points[-2], spline(0))
+    assert xp.allclose(points[-1], spline(spline.M - 1))
 
 
 def test_protected_spline_attributes(spline_curve, coeff_gen, basis_function):
@@ -1111,132 +1190,270 @@ def test_protected_spline_attributes(spline_curve, coeff_gen, basis_function):
         spline.pad = spline.pad + 1
 
 
-def test_single_value_input_call(initialized_spline_curve, single_value_t):
-    spline = initialized_spline_curve
-
+def test_single_value_input_call(initialised_spline_curve, single_value_t, xp, dtype):
+    spline = initialised_spline_curve
     output = spline(single_value_t)
-    assert isinstance(output, np.ndarray)
-    assert np.issubdtype(output.dtype, np.floating)
-    assert np.all(output.shape == (spline.ndim,))
+    assert array_api_compat.is_array_api_obj(output)
+    assert array_api_compat.array_namespace(output) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert output.dtype == expected_dtype
+    assert output.shape == (spline.ndim,)
 
 
-def test_array_input_call(initialized_spline_curve, array_t):
-    spline = initialized_spline_curve
+def test_array_input_call(initialised_spline_curve, array_t, xp, dtype):
+    spline = initialised_spline_curve
     output = spline(array_t)
-    assert isinstance(output, np.ndarray)
-    assert np.issubdtype(output.dtype, np.floating)
-    assert np.all(output.shape == (len(array_t), spline.ndim))
+    assert array_api_compat.is_array_api_obj(output)
+    assert array_api_compat.array_namespace(output) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert output.dtype == expected_dtype
+    assert output.shape == (len(array_t), spline.ndim)
 
 
-def test_single_value_input_arc_length(initialized_spline_curve, single_value_t):
-    spline = initialized_spline_curve
+def test_single_value_input_arc_length(initialised_spline_curve, single_value_t, xp, dtype):
+    spline = initialised_spline_curve
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
 
     arc_length = spline.arc_length(single_value_t)
-    assert isinstance(arc_length, float)
+    assert xp.isdtype(arc_length.dtype, expected_dtype)
 
     output = spline.arc_length_to_parameter(arc_length)
-    assert isinstance(output, float)
+    assert xp.isdtype(output.dtype, expected_dtype)
 
 
-def test_array_input_arc_length(initialized_spline_curve, array_t):
-    spline = initialized_spline_curve
+def test_array_input_arc_length(initialised_spline_curve, array_t, dtype, xp):
+    spline = initialised_spline_curve
 
     arc_lengths = spline.arc_length(array_t)
-    assert isinstance(arc_lengths, np.ndarray)
-    assert np.issubdtype(arc_lengths.dtype, np.floating)
-    assert np.all(arc_lengths.shape == (len(array_t),))
+    assert array_api_compat.is_array_api_obj(arc_lengths)
+    assert array_api_compat.array_namespace(arc_lengths) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert arc_lengths.dtype == expected_dtype
+    assert arc_lengths.shape == (len(array_t),)
 
     output = spline.arc_length_to_parameter(arc_lengths, atol=0.1)
-    assert isinstance(output, np.ndarray)
-    assert np.issubdtype(output.dtype, np.floating)
-    assert np.all(arc_lengths.shape == (len(array_t),))
+    assert array_api_compat.is_array_api_obj(output)
+    assert array_api_compat.array_namespace(output) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert output.dtype == expected_dtype
+    assert output.shape == (len(array_t),)
 
 
-def test_single_value_input_curvature(initialized_twice_differentiable_spline_curve, single_value_t):
-    spline = initialized_twice_differentiable_spline_curve
+def test_single_value_input_curvature(initialised_twice_differentiable_spline_curve, single_value_t, xp, dtype):
+    spline = initialised_twice_differentiable_spline_curve
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
 
     curvature = spline.curvature(single_value_t)
-    assert isinstance(curvature, float)
+    assert isinstance(curvature, expected_dtype)
 
 
-def test_array_input_curvature(initialized_twice_differentiable_spline_curve, array_t):
-    spline = initialized_twice_differentiable_spline_curve
+def test_array_input_curvature(initialised_twice_differentiable_spline_curve, array_t, xp, dtype):
+    spline = initialised_twice_differentiable_spline_curve
 
     curvatures = spline.curvature(array_t)
-    assert isinstance(curvatures, np.ndarray)
-    assert np.issubdtype(curvatures.dtype, np.floating)
-    assert np.all(curvatures.shape == (len(array_t),))
+    assert array_api_compat.is_array_api_obj(curvatures)
+    assert array_api_compat.array_namespace(curvatures) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert curvatures.dtype == expected_dtype
+    assert curvatures.shape == (len(array_t),)
 
 
-def test_single_value_input_dtheta(initialized_twice_differentiable_2D_spline_curve, single_value_t):
-    spline = initialized_twice_differentiable_2D_spline_curve
+def test_single_value_input_dtheta(initialised_twice_differentiable_2D_spline_curve, single_value_t, xp, dtype):
+    spline = initialised_twice_differentiable_2D_spline_curve
     dtheta = spline.dtheta(single_value_t)
-    assert isinstance(dtheta, float)
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert isinstance(dtheta, expected_dtype)
 
 
-def test_array_input_dtheta(initialized_twice_differentiable_2D_spline_curve, array_t):
-    spline = initialized_twice_differentiable_2D_spline_curve
+def test_array_input_dtheta(initialised_twice_differentiable_2D_spline_curve, array_t, xp, dtype):
+    spline = initialised_twice_differentiable_2D_spline_curve
 
     dthetas = spline.dtheta(array_t)
-    assert isinstance(dthetas, np.ndarray)
-    assert np.issubdtype(dthetas.dtype, np.floating)
-    assert np.all(dthetas.shape == (len(array_t),))
+    assert array_api_compat.is_array_api_obj(dthetas)
+    assert array_api_compat.array_namespace(dthetas) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert dthetas.dtype == expected_dtype
+    assert dthetas.shape == (len(array_t),)
 
 
-def test_single_value_input_normal_2D(initialized_2D_spline_curve, single_value_t):
-    spline = initialized_2D_spline_curve
+def test_single_value_input_normal_2D(initialised_2D_spline_curve, single_value_t, xp, dtype):
+    spline = initialised_2D_spline_curve
     output = spline.normal(single_value_t, shift=1e-10)
-    assert isinstance(output, np.ndarray)
-    assert np.issubdtype(output.dtype, np.floating)
-    assert np.all(output.shape == (2,))
+    assert array_api_compat.is_array_api_obj(output)
+    assert array_api_compat.array_namespace(output) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert output.dtype == expected_dtype
+    assert output.shape == (2,)
 
 
-def test_array_input_normal_2D(initialized_2D_spline_curve, array_t):
-    spline = initialized_2D_spline_curve
+def test_array_input_normal_2D(initialised_2D_spline_curve, array_t, xp, dtype):
+    spline = initialised_2D_spline_curve
     output = spline.normal(array_t, shift=1e-10)
-    assert isinstance(output, np.ndarray)
-    assert np.issubdtype(output.dtype, np.floating)
-    assert np.all(output.shape == (len(array_t), 2))
+    assert array_api_compat.is_array_api_obj(output)
+    assert array_api_compat.array_namespace(output) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert output.dtype == expected_dtype
+    assert output.shape == (len(array_t), 2)
 
 
-def test_single_value_input_normal_3D(initialized_twice_differentiable_3D_spline_curve, single_value_t, frame):
-    spline = initialized_twice_differentiable_3D_spline_curve
+def test_single_value_input_normal_3D(
+    initialised_twice_differentiable_3D_spline_curve, single_value_t, frame, xp, dtype
+):
+    spline = initialised_twice_differentiable_3D_spline_curve
     output = spline.normal(single_value_t, frame=frame, shift=1e-10)
-    assert isinstance(output, np.ndarray)
-    assert np.issubdtype(output.dtype, np.floating)
-    assert np.all(output.shape == (2, 3))
+    assert array_api_compat.is_array_api_obj(output)
+    assert array_api_compat.array_namespace(output) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert output.dtype == expected_dtype
+    assert output.shape == (2, 3)
 
 
-def test_array_input_normal_3D(initialized_twice_differentiable_3D_spline_curve, array_t, frame):
-    spline = initialized_twice_differentiable_3D_spline_curve
+def test_array_input_normal_3D(initialised_twice_differentiable_3D_spline_curve, array_t, frame, xp, dtype):
+    spline = initialised_twice_differentiable_3D_spline_curve
     output = spline.normal(array_t, frame=frame, shift=1e-10)
-    assert isinstance(output, np.ndarray)
-    assert np.issubdtype(output.dtype, np.floating)
-    assert np.all(output.shape == (len(array_t), 2, 3))
+    assert array_api_compat.is_array_api_obj(output)
+    assert array_api_compat.array_namespace(output) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert output.dtype == expected_dtype
+    assert output.shape == (len(array_t), 2, 3)
 
 
-def test_single_value_input_moving_frame_3D(initialized_twice_differentiable_3D_spline_curve, single_value_t, frame):
-    spline = initialized_twice_differentiable_3D_spline_curve
+def test_single_value_input_moving_frame_3D(
+    initialised_twice_differentiable_3D_spline_curve, single_value_t, frame, xp, dtype
+):
+    spline = initialised_twice_differentiable_3D_spline_curve
     output = spline.moving_frame(single_value_t, method=frame, shift=1e-10)
-    assert isinstance(output, np.ndarray)
-    assert np.issubdtype(output.dtype, np.floating)
-    assert np.all(output.shape == (3, 3))
+    assert array_api_compat.is_array_api_obj(output)
+    assert array_api_compat.array_namespace(output) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert output.dtype == expected_dtype
+    assert output.shape == (3, 3)
 
 
-def test_array_input_moving_frame_3D(initialized_twice_differentiable_3D_spline_curve, array_t, frame):
-    spline = initialized_twice_differentiable_3D_spline_curve
+def test_array_input_moving_frame_3D(initialised_twice_differentiable_3D_spline_curve, array_t, frame, xp, dtype):
+    spline = initialised_twice_differentiable_3D_spline_curve
     output = spline.moving_frame(array_t, method=frame, shift=1e-10)
-    assert isinstance(output, np.ndarray)
-    assert np.issubdtype(output.dtype, np.floating)
-    assert np.all(output.shape == (len(array_t), 3, 3))
+    assert array_api_compat.is_array_api_obj(output)
+    assert array_api_compat.array_namespace(output) is xp
+    expected_dtype = xp.asarray(1.0).dtype if dtype is int else dtype
+    assert output.dtype == expected_dtype
+    assert output.shape == (len(array_t), 3, 3)
 
 
-def test_ndim(initialized_spline_curve):
-    spline = initialized_spline_curve
-    expected = spline(np.array([0.5])).shape[-1]
+def test_ndim(initialised_spline_curve, as_xp_array):
+    spline = initialised_spline_curve
+    expected = spline(as_xp_array([0.5])).shape[-1]
     assert spline.ndim == expected
 
 
+def _basis_matrix_reference(spline, control_points_np):
+    """Construct a numpy spline with the same configuration for reference."""
+    return splinebox.spline_curves.Spline(
+        spline.M, spline.basis_function, closed=spline.closed, control_points=control_points_np
+    )
+
+
+def test_basis_matrix_matches_numpy_reference(spline_curve, as_xp_array, xp, is_hermite_spline, request, coeff_gen):
+    spline = spline_curve
+
+    if is_hermite_spline(spline):
+        # Currently not implemented for Hermite spline curves
+        request.node.add_marker(pytest.mark.xfail)
+
+    M = spline.M
+    support = spline.basis_function.support
+    closed = spline.closed
+    rng = np.random.default_rng(seed=1234)
+    n = M if closed else M + 2 * (math.ceil(support / 2) - 1)
+    control_points_np = rng.random((n, 3)) * 100
+
+    spline.control_points = as_xp_array(control_points_np)
+    t = as_xp_array(np.linspace(0, M - 1, 20))
+
+    bm = spline.basis_matrix(t)
+    result = bm @ spline.control_points
+
+    reference = _basis_matrix_reference(spline, control_points_np)
+    expected = reference(np.asarray(t))
+
+    assert np.allclose(np.asarray(result), expected)
+
+
+def test_basis_matrix_single_value_t(non_hermite_spline_curve, as_xp_array, xp):
+    spline = non_hermite_spline_curve
+    M = spline.M
+    support = spline.basis_function.support
+    closed = spline.closed
+    rng = np.random.default_rng(seed=1234)
+    n = M if closed else M + 2 * (math.ceil(support / 2) - 1)
+    control_points_np = rng.random((n, 3)) * 100
+
+    spline.control_points = as_xp_array(control_points_np)
+
+    bm = spline.basis_matrix(as_xp_array(np.array([1.5])))
+    assert bm.shape == (1, n)
+
+    result = bm @ spline.control_points
+    reference = _basis_matrix_reference(spline, control_points_np)
+    assert np.allclose(np.asarray(result[0]), reference(np.array(1.5))[0])
+
+
+def test_basis_matrix_out_of_range_t(non_hermite_spline_curve, as_xp_array, xp, closed):
+    """
+    The rows of the basis matrix correspond to the input t values, even if
+    some t values are outside the valid range of an open spline, and rows
+    without any valid control point are all zero.
+
+    The numpy path currently re-packs the rows in this case, so it is
+    skipped for open splines.
+    """
+    if not closed and xp.__name__.endswith("numpy"):
+        pytest.skip("The numpy path re-packs the rows for out-of-range t values.")
+    spline = non_hermite_spline_curve
+    M = spline.M
+    support = spline.basis_function.support
+    rng = np.random.default_rng(seed=1234)
+    n = M if closed else M + 2 * (math.ceil(support / 2) - 1)
+    control_points_np = rng.random((n, 3)) * 100
+
+    spline.control_points = as_xp_array(control_points_np)
+    # t values before and after the valid range of an open spline
+    t = np.linspace(-2, M + 2, 20)
+
+    bm = spline.basis_matrix(as_xp_array(t))
+    if scipy.sparse.issparse(bm):
+        bm = bm.toarray()
+    assert bm.shape == (len(t), n)
+
+    for i, t_value in enumerate(t):
+        row = np.asarray(bm[i])
+        single_row = spline.basis_matrix(as_xp_array(np.array([t_value])))
+        if scipy.sparse.issparse(single_row):
+            single_row = single_row.toarray()
+        assert np.allclose(row, np.asarray(single_row)[0], equal_nan=True)
+
+
+def test_basis_matrix_derivative_argument(non_hermite_spline_curve, as_xp_array, xp, derivative):
+    spline = non_hermite_spline_curve
+    M = spline.M
+    support = spline.basis_function.support
+    closed = spline.closed
+    rng = np.random.default_rng(seed=1234)
+    n = M if closed else M + 2 * (math.ceil(support / 2) - 1)
+    control_points_np = rng.random((n, 3)) * 100
+
+    spline.control_points = as_xp_array(control_points_np)
+    t = as_xp_array(np.linspace(0, M - 1, 20))
+
+    bm = spline.basis_matrix(t, derivative=derivative)
+    result = bm @ spline.control_points
+
+    reference = _basis_matrix_reference(spline, control_points_np)
+    expected = reference(np.asarray(t), derivative=derivative)
+
+    assert np.allclose(np.asarray(result), expected, equal_nan=True)
+
+
 def _test_t(spline, n=20):
     return np.linspace(0, spline.M, n, endpoint=False) if spline.closed else np.linspace(0, spline.M - 1, n)
 
@@ -1247,16 +1464,16 @@ def _basis_matrix_control_point_component(spline, t, derivative=0):
     return bm[0] if isinstance(bm, tuple) else bm
 
 
-def test_derivative_wrt_control_points_matches_basis_matrix(initialized_spline_curve):
-    spline = initialized_spline_curve
+def test_derivative_wrt_control_points_matches_basis_matrix(initialised_spline_curve, xp):
+    spline = initialised_spline_curve
     t = _test_t(spline)
     expected = _basis_matrix_control_point_component(spline, t).toarray()
     result = spline.derivative_wrt_control_points(t).toarray()
-    assert np.allclose(result, expected, equal_nan=True)
+    assert xp.allclose(result, expected, equal_nan=True)
 
 
-def test_derivative_wrt_control_points_finite_difference(initialized_spline_curve):
-    spline = initialized_spline_curve
+def test_derivative_wrt_control_points_finite_difference(initialised_spline_curve):
+    spline = initialised_spline_curve
     t = _test_t(spline)
 
     delta = np.random.randn(*spline.control_points.shape) * 1e-7
@@ -1269,16 +1486,16 @@ def test_derivative_wrt_control_points_finite_difference(initialized_spline_curv
     assert np.allclose(predicted, actual, atol=1e-5)
 
 
-def test_derivative_of_norm_squared_wrt_control_points_shape(initialized_spline_curve):
-    spline = initialized_spline_curve
+def test_derivative_of_norm_squared_wrt_control_points_shape(initialised_spline_curve):
+    spline = initialised_spline_curve
     n_control_points = spline.M if spline.closed else spline.M + 2 * spline.pad
     t = _test_t(spline)
     result = spline.derivative_of_norm_squared_wrt_control_points(t)
     assert result.shape == (len(t), n_control_points, spline.ndim)
 
 
-def test_derivative_of_norm_squared_wrt_control_points_finite_difference(initialized_spline_curve):
-    spline = initialized_spline_curve
+def test_derivative_of_norm_squared_wrt_control_points_finite_difference(initialised_spline_curve):
+    spline = initialised_spline_curve
     t = _test_t(spline)
 
     delta = np.random.randn(*spline.control_points.shape) * 1e-7
@@ -1295,8 +1512,8 @@ def total_norm_squared(s):
     assert np.isclose(predicted, actual, atol=1e-5)
 
 
-def test_derivative_wrt_control_points_derivative_argument(initialized_spline_curve):
-    spline = initialized_spline_curve
+def test_derivative_wrt_control_points_derivative_argument(initialised_spline_curve):
+    spline = initialised_spline_curve
     t = _test_t(spline)
     for derivative in [0, 1, 2]:
         expected = _basis_matrix_control_point_component(spline, t, derivative=derivative).toarray()
@@ -1330,13 +1547,13 @@ def test_derivative_wrt_tangents_finite_difference(hermite_spline_curve, coeff_g
     assert np.allclose(predicted, actual, atol=1e-5)
 
 
-def test_h1_finite_difference(initialized_spline_curve, is_hermite_spline, request):
+def test_h1_finite_difference(initialised_spline_curve, is_hermite_spline, request):
     # h1 appears in the derivative of the norm of the
     # first derivative of the spline to the fourth power.
     # For details, see the Active contour model sections of
     # in the Theory part of the documentation.
 
-    spline = initialized_spline_curve
+    spline = initialised_spline_curve
 
     if is_hermite_spline(spline):
         # Currently not implemented for Hermite spline curves
@@ -1362,13 +1579,13 @@ def test_h1_finite_difference(initialized_spline_curve, is_hermite_spline, reque
     assert np.allclose(predicted, actual, atol=1e-3)
 
 
-def test_h2_finite_difference(initialized_spline_curve, is_hermite_spline, request):
+def test_h2_finite_difference(initialised_spline_curve, is_hermite_spline, request):
     # h2 appears in the derivative of the squared norm of the
     # first derivative of the spline. For details, see the
     # Active contour model sections of in the Theory part of the
     # documentation.
 
-    spline = initialized_spline_curve
+    spline = initialised_spline_curve
 
     if is_hermite_spline(spline):
         # Currently not implemented for Hermite spline curves
@@ -1420,9 +1637,9 @@ def _numerical_curvature_energy(spline, n=10000):
 
 
 def test_derivative_of_curvilinear_reparametrisation_energy_wrt_control_points_shape(
-    initialized_spline_curve, is_hermite_spline, request
+    initialised_spline_curve, is_hermite_spline, request
 ):
-    spline = initialized_spline_curve
+    spline = initialised_spline_curve
     if is_hermite_spline(spline):
         # Currently not implemented for Hermite spline curves
         request.node.add_marker(pytest.mark.xfail)
@@ -1490,9 +1707,9 @@ def test_derivative_of_curvilinear_reparametrisation_energy_zero_for_uniform_spe
 
 
 def test_derivative_of_curvilinear_reparametrisation_energy_wrt_control_points_finite_difference(
-    initialized_spline_curve, is_hermite_spline, request
+    initialised_spline_curve, is_hermite_spline, request
 ):
-    spline = initialized_spline_curve
+    spline = initialised_spline_curve
     if is_hermite_spline(spline):
         # Currently not implemented for Hermite spline curves
         request.node.add_marker(pytest.mark.xfail)
@@ -1512,9 +1729,9 @@ def test_derivative_of_curvilinear_reparametrisation_energy_wrt_control_points_f
 
 
 def test_derivative_of_curvilinear_reparametrisation_energy_wrt_control_points_default_c(
-    initialized_spline_curve, is_hermite_spline, request
+    initialised_spline_curve, is_hermite_spline, request
 ):
-    spline = initialized_spline_curve
+    spline = initialised_spline_curve
     if is_hermite_spline(spline):
         # Currently not implemented for Hermite spline curves
         request.node.add_marker(pytest.mark.xfail)
@@ -1540,16 +1757,8 @@ def test_derivative_of_curvilinear_reparametrisation_energy_wrt_control_points_d
     assert np.isclose(predicted, actual, rtol=1e-3, atol=1e-9)
 
 
-# def test_derivative_of_curvilinear_reparametrisation_energy_not_implemented_for_hermite_spline(
-#     initialized_hermite_spline_curve,
-# ):
-#     spline = initialized_hermite_spline_curve
-#     with pytest.raises(RuntimeError):
-#         spline.derivative_of_curvilinear_reparametrisation_energy_wrt_control_points(1.0)
-
-
-def test_derivative_of_curvature_energy_wrt_control_points_shape(initialized_spline_curve, is_hermite_spline, request):
-    spline = initialized_spline_curve
+def test_derivative_of_curvature_energy_wrt_control_points_shape(initialised_spline_curve, is_hermite_spline, request):
+    spline = initialised_spline_curve
     if is_hermite_spline(spline):
         # Currently not implemented for Hermite spline curves
         request.node.add_marker(pytest.mark.xfail)
@@ -1599,9 +1808,9 @@ def test_derivative_of_curvature_energy_wrt_control_points_radial_for_circle(M,
 
 
 def test_derivative_of_curvature_energy_wrt_control_points_finite_difference(
-    initialized_spline_curve, is_hermite_spline, request
+    initialised_spline_curve, is_hermite_spline, request
 ):
-    spline = initialized_spline_curve
+    spline = initialised_spline_curve
     if is_hermite_spline(spline):
         # Currently not implemented for Hermite spline curves
         request.node.add_marker(pytest.mark.xfail)
diff --git a/tox.ini b/tox.ini
index 7d2958f..50dd293 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,5 +1,5 @@
 [tox]
-envlist = py{39,310,311,312,313}-{linux,macos,windows}
+envlist = py{39,310,311,312,313,314}-{linux,macos,windows}
 isolated_build=true
 
 [gh-actions]
@@ -9,6 +9,7 @@ python =
     3.11: py311
     3.12: py312
     3.13: py313
+    3.14: py314
 
 [gh-actions:env]
 PLATFORM =
@@ -24,8 +25,16 @@ platform =
 passenv =
     CI
     GITHUB_ACTIONS
+    USERNAME
+    USER
+    LOGNAME
+    PATH
+    INCLUDE
+    LIB
+    LIBPATH
+    TORCHINDUCTOR_CACHE_DIR
 extras =
     test,examples
 commands =
-    pytest -v --color=yes --doctest-modules src
-    pytest -v --color=yes --cov=splinebox --cov-report=xml
+    pytest -vx --color=yes --doctest-modules src
+    pytest -vx --color=yes --cov=splinebox --cov-report=xml