diff --git a/.github/gpu_ci/README.md b/.github/gpu_ci/README.md
index 7466aa05e..91857ded4 100644
--- a/.github/gpu_ci/README.md
+++ b/.github/gpu_ci/README.md
@@ -22,7 +22,7 @@ allocation per workflow run, and every CI job of that run executes in it:
 | `build_madspace` | cluster | build madspace with `ENABLE_CUDA`/`ENABLE_HIP` for the GPU of the node (`build_madspace_gpu.sh`). The build is kept in `$GLOBALSCRATCH/mg7-gpu-ci/cache` and redone only when the madspace sources, the modules or the GPU architecture change |
 | `pp_ttx` | cluster | `generate p p > t t~`, `output mg7`, `device = ["cuda"]` or `["hip"]` in `run_card.toml`, `bin/generate_events -f`, check the cross section in `info.json` (`pp_ttx_mg7.sh`). The logs and cards are uploaded as an artifact, and the cross section is shown in the run summary |
 | `pp_jj` | cluster | the same with `generate p p > j j` (`PROCESS`/`TAG` of `pp_ttx_mg7.sh`): its jet cuts act on several objects, which `p p > t t~` has none of |
-| `madspace_ops` | cluster | madspace GPU runtime operations against the CPU runtime on the same inputs (`madspace_ops_gpu.sh`, `madspace_checks.py`; torch puts the inputs on the GPU): a histogram with underflow and overflow bins, an elementwise kernel on a rank-4 tensor |
+| `madspace_tests` | cluster | the madspace unit tests, `madspace/tests`, with `MADSPACE_TEST_DEVICE=cuda` or `hip` (`madspace_tests_gpu.sh`): the tests put their torch tensors on the GPU, so madspace runs them on its GPU runtime. Host-side tests (marker `cpu_only`) are skipped. pytest and scipy are installed in the cache when the modules lack them, and so is a CUDA torch wheel when the PyTorch module is a CPU build (manneback). The log and `junit.xml` are uploaded as an artifact |
 | `interference` | cluster | squared split orders on the GPU backend of madmatrix (`interference_gpu.sh`, `interference_checks.py`): standalone `u u~ > u u~` constraints, the helicity choice of an interference matrix element through umami (exact and summed), and `p p > u u~ QCD^2==2` with mg7 against its CPU reference. Input `interference`, on by default and on push |
 | `stop_runner` | cluster | clean up, then create the stop file: `runner_batch.sh` stops the runner and the allocation ends |
 | `watchdog` | vmadgraph | from the start of the runner to the end of the allocation: if the allocation ends before `stop_runner` completed, cancel the run (`ceci-gate watch-runner`, `remote_watch_runner.sh`) |
diff --git a/.github/gpu_ci/madspace_checks.py b/.github/gpu_ci/madspace_checks.py
deleted file mode 100644
index c2a557b65..000000000
--- a/.github/gpu_ci/madspace_checks.py
+++ /dev/null
@@ -1,109 +0,0 @@
-#!/usr/bin/env python3
-"""madspace GPU runtime operations against the CPU runtime, on the same inputs.
-
-Run by madspace_ops_gpu.sh (job madspace_ops of gpu_runner_ci.yml) in its run
-directory. The inputs go to the GPU as torch tensors (dlpack), so that madspace runs
-the function on its GPU runtime; the same function on numpy inputs runs on the CPU.
-
-1. histogram: a weighted histogram with underflow and overflow bins, its values and
-   squared values (the GPU reduction used to write its n_bins + 2 keys into a scratch
-   buffer of n_bins, and summed the squared weights of the wrong events in each bin);
-2. rank4_add: an elementwise kernel on a rank-4 tensor, batch x 2 x 3 x 4 (the GPU
-   launcher computed only the first slice of the last dimension).
-
-Without a GPU-enabled torch among the modules (manneback: the default PyTorch module is
-the CPU build) the checks are reported as skipped. Writes summary.txt (key=value lines)
-in the current directory; exit status 1 if a check fails.
-"""
-
-import os
-import sys
-import time
-
-import numpy as np
-
-import madspace as ms
-
-results = {}
-
-
-def check(name, ok, detail):
-    print('%s %s: %s' % ('PASS' if ok else 'FAIL', name, detail), flush=True)
-    results[name] = (ok, detail)
-
-
-def to_numpy(tensor, torch):
-    return torch.from_dlpack(tensor).cpu().numpy()
-
-
-def check_histogram(torch, device, n=100000, n_bins=40):
-    out_type = ms.Type(ms.DataType.float, ms.BatchSize.one, [n_bins + 2])
-    fb = ms.FunctionBuilder(
-        ms.NamedTypes([(name, ms.batch_float) for name in ('x', 'w', 'lo', 'hi')]),
-        ms.NamedTypes([('values', out_type), ('square_values', out_type)]),
-    )
-    values, square_values = fb.histogram(*(fb.input(i) for i in range(4)), n_bins)
-    fb.output(0, values)
-    fb.output(1, square_values)
-    function = fb.function()
-    rng = np.random.default_rng(17)
-    inputs = [rng.uniform(-0.5, 1.5, n), rng.normal(size=n), np.zeros(n), np.ones(n)]
-    cpu = [np.asarray(ms.Tensor.numpy(t))
-           for t in ms.FunctionRuntime(function, ms.default_context()).call(inputs)]
-    gpu = [to_numpy(t, torch) for t in ms.FunctionRuntime(function).call(
-        [torch.from_numpy(a).to(device) for a in inputs])]
-    worst = [float(np.max(np.abs(g - c) / np.maximum(np.abs(c), 1e-300)))
-             for g, c in zip(gpu, cpu)]
-    under, over = cpu[0][0][0], cpu[0][0][-1]
-    check('histogram', max(worst) < 1e-9 and under != 0 and over != 0,
-          '%d bins + underflow/overflow, max rel diff GPU vs CPU %.1e (values), %.1e '
-          '(squared values); underflow %.4g, overflow %.4g'
-          % (n_bins, worst[0], worst[1], under, over))
-
-
-def check_rank4_add(torch, device, n=4096):
-    rank4 = ms.Type(ms.DataType.float, ms.batch_size, [2, 3, 4])
-    fb = ms.FunctionBuilder(ms.NamedTypes([('x', rank4)]), ms.NamedTypes([('y', rank4)]))
-    x = fb.input(0)
-    fb.output(0, fb.add(x, x))
-    a = np.random.default_rng(23).normal(size=(n, 2, 3, 4))
-    (y,) = ms.FunctionRuntime(fb.function()).call([torch.from_numpy(a).to(device)])
-    y = to_numpy(y, torch)
-    wrong = int(np.sum(y != 2 * a))
-    check('rank4_add', wrong == 0, '%d of %d elements of x + x differ from 2 x'
-          % (wrong, a.size))
-
-
-def main():
-    start = time.time()
-    try:
-        import torch
-    except ImportError as error:
-        check('torch', False, 'cannot import torch (%s): no GPU inputs' % error)
-        torch = None
-    skipped = None
-    if torch is not None and not torch.cuda.is_available():
-        skipped = 'torch %s of the modules has no GPU support' % torch.__version__
-        print('SKIP all checks: %s' % skipped, flush=True)
-        torch = None
-    if torch is not None:
-        for name, step in (('histogram', check_histogram), ('rank4_add', check_rank4_add)):
-            try:
-                step(torch, 'cuda')  # also the ROCm devices in torch
-            except Exception as error:  # one broken step must not hide the others
-                check(name, False, 'error: %s' % error)
-    failures = [name for name, (ok, _detail) in results.items() if not ok]
-    with open('summary.txt', 'w') as f:
-        for name, (ok, detail) in results.items():
-            f.write('%s=%s (%s)\n' % (name, 'ok' if ok else 'FAILED', detail))
-        if skipped:
-            f.write('checks=skipped (%s)\n' % skipped)
-        f.write('node=%s\ngpu=%s\nwalltime=%ds\n' % (
-            os.uname().nodename, os.environ.get('GPU_NAME', 'none'), time.time() - start))
-    print('\n%d check(s) failed: %s' % (len(failures), ', '.join(failures)) if failures
-          else '\nskipped: %s' % skipped if skipped else '\nall checks passed')
-    return 1 if failures else 0
-
-
-if __name__ == '__main__':
-    sys.exit(main())
diff --git a/.github/gpu_ci/madspace_ops_gpu.sh b/.github/gpu_ci/madspace_ops_gpu.sh
deleted file mode 100755
index 3f30a6e6a..000000000
--- a/.github/gpu_ci/madspace_ops_gpu.sh
+++ /dev/null
@@ -1,33 +0,0 @@
-#!/bin/bash -l
-# Runs on the self-hosted runner (GPU node): job madspace_ops of gpu_runner_ci.yml.
-# madspace GPU runtime operations against the CPU runtime: see madspace_checks.py.
-# Environment: as pp_ttx_mg7.sh (BACKEND, MODULES, MADSPACE_PREFIX, WORKDIR).
-# torch (from the PyTorch module) puts the inputs on the GPU. This runs the python3 of
-# the modules, not the venv of build_madspace (same interpreter, which madspace was built
-# with): EasyBuild finds some packages of the modules (e.g. typing_extensions of
-# Python-bundle-PyPI, which torch imports) through EBPYTHONPREFIXES, which only the
-# python3 of the modules reads (its sitecustomize), not a venv made from it.
-set -eo pipefail
-HERE=$(cd "$(dirname "$0")" && pwd)
-section() { echo; echo "=================== $* ($(date +%T))"; }
-
-section "Environment"
-if [ -n "$MODULES" ]; then source "$HERE/load_modules.sh"; fi
-python3 --version
-case $BACKEND in
-    cuda) GPU_NAME=$(nvidia-smi --query-gpu=name --format=csv,noheader | sed -n 1p) ;;
-    hip)  GPU_NAME=$(rocminfo 2> /dev/null | awk -F': *' '/Marketing Name/ {name=$2} /Device Type: *GPU/ {print name; exit}' || true) ;;
-    *)    echo "::error::the madspace checks need a GPU backend, not '$BACKEND'"; exit 2 ;;
-esac
-export GPU_NAME
-export PYTHONPATH=$MADSPACE_PREFIX${PYTHONPATH:+:$PYTHONPATH}
-
-rm -rf "$WORKDIR"
-mkdir -p "$WORKDIR"
-cd "$WORKDIR"
-STATUS=0
-python3 "$HERE/madspace_checks.py" 2>&1 | tee checks.log || STATUS=1
-
-section "Summary"
-cat summary.txt 2> /dev/null || echo "no summary.txt"
-exit $STATUS
diff --git a/.github/gpu_ci/madspace_tests_gpu.sh b/.github/gpu_ci/madspace_tests_gpu.sh
new file mode 100755
index 000000000..db942911a
--- /dev/null
+++ b/.github/gpu_ci/madspace_tests_gpu.sh
@@ -0,0 +1,65 @@
+#!/bin/bash -l
+# Runs on the self-hosted runner (GPU node): job madspace_tests of gpu_runner_ci.yml.
+# The madspace unit tests (madspace/tests) on the GPU: MADSPACE_TEST_DEVICE=$BACKEND puts
+# their torch tensors on the GPU, so that madspace runs on its GPU runtime.
+# Environment: as pp_ttx_mg7.sh (BACKEND, MODULES, MADSPACE_PREFIX, WORKDIR).
+# This runs the python3 of the modules, not the venv of build_madspace (same interpreter,
+# which madspace was built with): EasyBuild finds some packages of the modules (e.g.
+# typing_extensions of Python-bundle-PyPI, which torch imports) through EBPYTHONPREFIXES,
+# which only the python3 of the modules reads (its sitecustomize), not a venv made from it.
+set -eo pipefail
+HERE=$(cd "$(dirname "$0")" && pwd)
+REPO=$(cd "$HERE/../.." && pwd)
+section() { echo; echo "=================== $* ($(date +%T))"; }
+START=$(date +%s)
+
+section "Environment"
+if [ -n "$MODULES" ]; then source "$HERE/load_modules.sh"; fi
+python3 --version
+case $BACKEND in
+    cuda) GPU_NAME=$(nvidia-smi --query-gpu=name --format=csv,noheader | sed -n 1p) ;;
+    hip)  GPU_NAME=$(rocminfo 2> /dev/null | awk -F': *' '/Marketing Name/ {name=$2} /Device Type: *GPU/ {print name; exit}' || true) ;;
+    *)    echo "::error::the madspace tests on the GPU need a GPU backend, not '$BACKEND'"; exit 2 ;;
+esac
+# test requirements the modules lack (pytest, scipy) are kept per python version, after
+# the packages of the modules in the path
+DEPS=${GLOBALSCRATCH:-$HOME}/mg7-gpu-ci/cache/test-deps-$(python3 -c 'import sys; print("%d.%d" % sys.version_info[:2])')
+export PYTHONPATH=$MADSPACE_PREFIX${PYTHONPATH:+:$PYTHONPATH}:$DEPS
+MISSING=$(python3 -c 'import importlib.util as u; print(" ".join(m for m in ("pytest", "scipy") if u.find_spec(m) is None))')
+if [ -n "$MISSING" ]; then
+    echo "Installing $MISSING into $DEPS"
+    python3 -m pip install --quiet --target="$DEPS" $MISSING
+fi
+# manneback has only a CPU build of PyTorch among its modules: then a CUDA wheel, kept in
+# the cache, goes before it in the path
+if [ "$BACKEND" = cuda ] && ! python3 -c 'import sys, torch; sys.exit(not torch.cuda.is_available())' 2> /dev/null; then
+    TORCH_WHEEL=${TORCH_WHEEL:-torch==2.8.0}
+    TORCH_INDEX=${TORCH_INDEX:-https://download.pytorch.org/whl/cu128}
+    TORCH_DIR=$DEPS-$(printf '%s %s' "$TORCH_WHEEL" "$TORCH_INDEX" | sha1sum | cut -c1-8)
+    if [ ! -e "$TORCH_DIR/.ready" ]; then
+        echo "The torch of the modules has no CUDA: installing $TORCH_WHEEL from $TORCH_INDEX into $TORCH_DIR"
+        rm -rf "$TORCH_DIR"
+        python3 -m pip install --quiet --target="$TORCH_DIR" --index-url "$TORCH_INDEX" "$TORCH_WHEEL"
+        touch "$TORCH_DIR/.ready"
+    fi
+    export PYTHONPATH=$TORCH_DIR:$PYTHONPATH
+fi
+python3 -c 'import torch; print("torch", torch.__version__, "GPU available:", torch.cuda.is_available())'
+
+rm -rf "$WORKDIR"
+mkdir -p "$WORKDIR"
+cd "$WORKDIR"
+section "madspace tests on $BACKEND ($GPU_NAME)"
+STATUS=0
+MADSPACE_TEST_DEVICE=$BACKEND python3 -m pytest "$REPO/madspace/tests" -ra -p no:cacheprovider \
+    --junitxml=junit.xml 2>&1 | tee tests.log || STATUS=1
+
+section "Summary"
+{
+    echo "tests=$(grep -E '^=+ .*(passed|failed|error|skipped).* =+$' tests.log | tail -1 | sed 's/^=* *//; s/ *=*$//')"
+    echo "node=$(uname -n)"
+    echo "gpu=${GPU_NAME:-none}"
+    echo "walltime=$(( $(date +%s) - START ))s"
+} > summary.txt
+cat summary.txt
+exit $STATUS
diff --git a/.github/workflows/gpu_runner_ci.yml b/.github/workflows/gpu_runner_ci.yml
index b705f2ba1..fb2a5490c 100644
--- a/.github/workflows/gpu_runner_ci.yml
+++ b/.github/workflows/gpu_runner_ci.yml
@@ -7,7 +7,7 @@
 #   build_madspace  cluster runner  madspace with the CUDA/HIP backend (cached on the cluster)
 #   pp_ttx          cluster runner  p p > t t~ with mg7 on the GPU
 #   pp_jj           cluster runner  p p > j j with mg7 on the GPU (cuts on several jets)
-#   madspace_ops    cluster runner  madspace GPU operations against the CPU runtime, .github/gpu_ci/madspace_checks.py
+#   madspace_tests  cluster runner  the madspace unit tests on the GPU, .github/gpu_ci/madspace_tests_gpu.sh
 #   interference    cluster runner  squared split orders (interference), .github/gpu_ci/interference_gpu.sh
 #   ...             cluster runner  more CI jobs: same needs, if and runs-on, and add them to
 #                                   stop_runner's needs
@@ -301,7 +301,8 @@ jobs:
           path: results/
           if-no-files-found: warn
 
-  madspace_ops:
+
+  madspace_tests:
     # madspace GPU runtime operations against the CPU runtime (histogram, rank-4 kernels)
     needs: [start_runner, build_madspace]
     if: needs.start_runner.outputs.attempt == github.run_attempt
@@ -313,42 +314,48 @@ jobs:
     env:
       BACKEND: ${{ inputs.backend }}
       MODULES: ${{ inputs.modules }}
-      # no VENV: the python3 of the modules (see madspace_ops_gpu.sh)
+      # no VENV: the python3 of the modules (see madspace_tests_gpu.sh)
       MADSPACE_PREFIX: ${{ needs.build_madspace.outputs.madspace_prefix }}
     steps:
       - uses: actions/checkout@v5
       - name: run directory
-        run: echo "WORKDIR=${GLOBALSCRATCH:-$HOME}/mg7-gpu-ci/runs/$RUN_TAG/madspace_ops" >> "$GITHUB_ENV"
-      - name: madspace operations on the GPU
-        run: .github/gpu_ci/madspace_ops_gpu.sh
+        run: echo "WORKDIR=${GLOBALSCRATCH:-$HOME}/mg7-gpu-ci/runs/$RUN_TAG/madspace_tests" >> "$GITHUB_ENV"
+      - name: madspace tests on the GPU
+        run: .github/gpu_ci/madspace_tests_gpu.sh
       - name: collect logs
         if: always()
         run: |
           mkdir -p results
-          (cd "$WORKDIR" && tar -cf - --ignore-failed-read summary.txt checks.log 2> /dev/null) | tar -xf - -C results || true
+          (cd "$WORKDIR" && tar -cf - --ignore-failed-read summary.txt tests.log junit.xml 2> /dev/null) | tar -xf - -C results || true
       - name: job summary
         if: always()
         run: |
           {
-            echo "### madspace operations on ${{ inputs.cluster }} ($BACKEND)"
+            echo "### madspace tests on ${{ inputs.cluster }} ($BACKEND)"
             if [ -f results/summary.txt ]; then
               echo "| | |"
               echo "|---|---|"
               sed 's/^\([^=]*\)=\(.*\)$/| \1 | \2 |/' results/summary.txt
             else
-              echo "No summary.txt: the job stopped early, see the log of the 'madspace operations on the GPU' step."
+              echo "No summary.txt: the job stopped early, see the log of the 'madspace tests on the GPU' step."
+            fi
+            if [ -f results/tests.log ]; then
+              echo
+              echo '```'
+              grep -E '^(FAILED|ERROR|XFAIL|SKIPPED) ' results/tests.log | head -100 || true
+              echo '```'
             fi
           } >> "$GITHUB_STEP_SUMMARY"
       - uses: actions/upload-artifact@v4
         if: always()
         with:
-          name: gpu-madspace_ops-${{ inputs.cluster }}-${{ inputs.backend }}
+          name: gpu-madspace_tests-${{ inputs.cluster }}-${{ inputs.backend }}
           path: results/
           if-no-files-found: warn
 
   stop_runner:
     # must list every job that runs on the self-hosted runner
-    needs: [start_runner, build_madspace, pp_ttx, pp_jj, madspace_ops, interference]
+    needs: [start_runner, build_madspace, pp_ttx, pp_jj, madspace_tests, interference]
     if: always() && needs.start_runner.result == 'success'
     # In an attempt that did not start the runner (a re-run of failed jobs, or of one job), a job
     # for the runner label would wait 24 hours, holding the concurrency group of the cluster.
diff --git a/.github/workflows/madspace.yml b/.github/workflows/madspace.yml
index 9044c0341..e499a6b5c 100644
--- a/.github/workflows/madspace.yml
+++ b/.github/workflows/madspace.yml
@@ -87,6 +87,8 @@ jobs:
               scipy
               --extra-index-url https://download.pytorch.org/whl/cpu
               torch
+            # the GPU runs of these tests are in gpu_runner_ci.yml (job madspace_tests)
+            CIBW_TEST_ENVIRONMENT: MADSPACE_TEST_DEVICE=cpu
             CIBW_TEST_COMMAND: pytest {project}/madspace/tests
 
       - uses: actions/upload-artifact@v4
diff --git a/madspace/src/gpu/runtime.cu b/madspace/src/gpu/runtime.cu
index 8e9a06a21..1b2624b72 100644
--- a/madspace/src/gpu/runtime.cu
+++ b/madspace/src/gpu/runtime.cu
@@ -40,8 +40,8 @@ void op_matrix_element(
     std::size_t output_count = locals[instruction.input_indices[2]].index_value();
     TensorVec contiguous_inputs(input_count);
     std::vector<UmamiInputKey> input_keys(input_count);
-    // one more output: the stream (UMAMI_OUT_GPU_STREAM), so that the matrix element runs
-    // in stream order with the rest of the function, not on the default stream
+    // one more output: the stream (UMAMI_OUT_GPU_STREAM), so that the matrix element
+    // runs in stream order with the rest of the function, not on the default stream
     std::vector<UmamiOutputKey> output_keys(output_count + 1);
     std::vector<void*> input_ptrs(input_count), output_ptrs(output_count + 1);
     for (std::size_t i = 0; i < input_count; ++i) {
@@ -1408,7 +1408,8 @@ void op_histogram(
     Tensor square_weights_tmp(
         DataType::dt_float, {padded_size}, device, AllocHint::temporary
     );
-    // the keys reduce_by_key writes: the n_bins bins, the underflow and the overflow bin
+    // the keys reduce_by_key writes: the n_bins bins, the underflow and the overflow
+    // bin
     Tensor reduce_tmp(DataType::dt_float, {n_bins + 2}, device, AllocHint::temporary);
 
     launch_kernel(
@@ -1583,7 +1584,8 @@ GpuRuntime::GpuRuntime(const Function& function_arg, ContextPtr context) :
     LastUseOfLocals last_use(function);
     InstructionDependencies dependencies(function);
 
-    std::size_t stream_count = 0, event_count = 0, backward_event_count = 0;
+    // at least the main stream, also for a function without instructions
+    std::size_t stream_count = 1, event_count = 0, backward_event_count = 0;
     for (auto& instr : function.instructions()) {
         if (instr.stream_index >= stream_count) {
             stream_count = instr.stream_index + 1;
diff --git a/madspace/src/kernels/kinematics.hpp b/madspace/src/kernels/kinematics.hpp
index b2ae2aba5..011e07dc6 100644
--- a/madspace/src/kernels/kinematics.hpp
+++ b/madspace/src/kernels/kinematics.hpp
@@ -355,8 +355,7 @@ KERNELSPEC FourMom<T> boost(FourMom<T> k, FourMom<T> p_boost, FVal<T> sign) {
 }
 
 template <typename T>
-KERNELSPEC FourMom<T>
-boost_light_cone(
+KERNELSPEC FourMom<T> boost_light_cone(
     FourMom<T> k, FVal<T> k_mass2, FourMom<T> p_boost, FVal<T> p_mass2, FVal<T> sign
 ) {
     // Same boost as `boost`, for a momentum k of known mass^2 k_mass2, written
@@ -389,12 +388,10 @@ boost_light_cone(
     auto k_plus_direct = k[0] + k_par;
     auto k_minus_direct = k[0] - k_par;
     auto forward = k_par >= 0.;
-    auto k_plus = where(
-        forward, k_plus_direct, trans2 / max(k_minus_direct, FVal<T>(EPS2))
-    );
-    auto k_minus = where(
-        forward, trans2 / max(k_plus_direct, FVal<T>(EPS2)), k_minus_direct
-    );
+    auto k_plus =
+        where(forward, k_plus_direct, trans2 / max(k_minus_direct, FVal<T>(EPS2)));
+    auto k_minus =
+        where(forward, trans2 / max(k_plus_direct, FVal<T>(EPS2)), k_minus_direct);
 
     auto k_plus_new = exp_y * k_plus;
     auto k_minus_new = exp_my * k_minus;
@@ -435,8 +432,7 @@ KERNELSPEC Pair<FourMom<T>, FourMom<T>> boost_two_body(
 }
 
 template <typename T>
-KERNELSPEC void
-boost_beam(
+KERNELSPEC void boost_beam(
     FIn<T, 2> q,
     FIn<T, 1> masses,
     FVal<T> x1,
@@ -476,12 +472,10 @@ boost_beam(
         auto trans2 = m2 + pt2;
         auto q_plus_direct = q0 + q3;
         auto q_minus_direct = q0 - q3;
-        auto q_plus = where(
-            q3 >= 0., q_plus_direct, trans2 / max(q_minus_direct, FVal<T>(EPS2))
-        );
-        auto q_minus = where(
-            q3 < 0., q_minus_direct, trans2 / max(q_plus_direct, FVal<T>(EPS2))
-        );
+        auto q_plus =
+            where(q3 >= 0., q_plus_direct, trans2 / max(q_minus_direct, FVal<T>(EPS2)));
+        auto q_minus =
+            where(q3 < 0., q_minus_direct, trans2 / max(q_plus_direct, FVal<T>(EPS2)));
         auto p_plus = e_plus * q_plus;
         auto p_minus = e_minus * q_minus;
         // Only for a momentum that is on its mass shell to rounding. The
@@ -491,8 +485,8 @@ boost_beam(
         // gained, and it keeps the linear boost.
         auto on_shell =
             fabs(q_plus_direct * q_minus_direct - trans2) <= 1e-13 * q0 * q0;
-        auto light_cone = (q_plus_direct >= 0.) & (q_minus_direct >= 0.) & on_shell &
-            (lin_0 < q0);
+        auto light_cone =
+            (q_plus_direct >= 0.) & (q_minus_direct >= 0.) & on_shell & (lin_0 < q0);
         p_out_i[0] = where(light_cone, 0.5 * (p_plus + p_minus), lin_0);
         p_out_i[1] = q1;
         p_out_i[2] = q2;
@@ -522,10 +516,10 @@ boost_beam(
             auto lin_0 = q0 * cosh_rap + q3 * sinh_rap;
             auto lin_3 = q3 * cosh_rap + q0 * sinh_rap;
             auto is_best = best_i == IVal<T>(static_cast<int>(i));
-            shift_0 = shift_0 +
-                where(is_best, FVal<T>(0.), FVal<T>(p_out[i][0]) - lin_0);
-            shift_3 = shift_3 +
-                where(is_best, FVal<T>(0.), FVal<T>(p_out[i][3]) - lin_3);
+            shift_0 =
+                shift_0 + where(is_best, FVal<T>(0.), FVal<T>(p_out[i][0]) - lin_0);
+            shift_3 =
+                shift_3 + where(is_best, FVal<T>(0.), FVal<T>(p_out[i][3]) - lin_3);
             best_lin_0 = where(is_best, lin_0, best_lin_0);
             best_lin_3 = where(is_best, lin_3, best_lin_3);
         }
@@ -562,10 +556,17 @@ KERNELSPEC Pair<FourMom<T>, FVal<T>> p1com_from_tabs_cs(
     FourMom<T> p1_com;
     auto e1_com = 0.5 * (m_tot + ed);
     p1_com[0] = max(e1_com, 0.);
-    p1_com[3] =
-        -(m1 * m1 + ma_2 + t_abs - 2. * p1_com[0] * pa_com[0]) / (2.0 * pa_com_mag);
-    auto pt2 = pp * pp - p1_com[3] * p1_com[3];
-    auto pt = sqrt(max(pt2, 0.));
+    // p1 along pa at |t| = t_lo, against it at t_lo + 4 |pa| pp. t_lo through
+    // E - p = m^2 / (E + p), and p1_z, p1_t from the distances to the ends of
+    // the range: pp^2 - p1_z^2 cancels near an end, and differently with FMA
+    auto plus_a = pa_com[0] + pa_com_mag;
+    auto plus_1 = p1_com[0] + pp;
+    auto t_lo = (plus_a - plus_1) * (m1 * m1 / plus_1 - ma_2 / plus_a);
+    auto pa_2 = 2. * pa_com_mag;
+    auto to_lo = (t_abs - t_lo) / pa_2;
+    auto to_hi = (t_lo + 2. * pa_2 * pp - t_abs) / pa_2;
+    p1_com[3] = 0.5 * (to_hi - to_lo);
+    auto pt = sqrt(max(to_lo * to_hi, 0.));
     p1_com[1] = pt * cos_phi;
     p1_com[2] = pt * sin_phi;
 
@@ -603,16 +604,14 @@ KERNELSPEC Quartuplet<FVal<T>, FVal<T>, FVal<T>, FVal<T>> phi_m1_m2_from_p1com(
 // Kernels
 
 template <typename T>
-KERNELSPEC void
-kernel_boost_beam(
+KERNELSPEC void kernel_boost_beam(
     FIn<T, 2> p1, FIn<T, 1> masses, FIn<T, 0> x1, FIn<T, 0> x2, FOut<T, 2> p_out
 ) {
     boost_beam<T>(p1, masses, x1, x2, 1.0, p_out);
 }
 
 template <typename T>
-KERNELSPEC void
-kernel_boost_beam_inverse(
+KERNELSPEC void kernel_boost_beam_inverse(
     FIn<T, 2> p1, FIn<T, 1> masses, FIn<T, 0> x1, FIn<T, 0> x2, FOut<T, 2> p_out
 ) {
     boost_beam<T>(p1, masses, x1, x2, -1.0, p_out);
diff --git a/madspace/tests/conftest.py b/madspace/tests/conftest.py
new file mode 100644
index 000000000..6dd5bbfff
--- /dev/null
+++ b/madspace/tests/conftest.py
@@ -0,0 +1,36 @@
+"""Runs the tests on the device set by MADSPACE_TEST_DEVICE (cpu, cuda or hip).
+
+Tests marked cpu_only (host-side code, e.g. LHE files) are skipped on a GPU.
+"""
+
+import pytest
+import torch
+from madspace_testing import BACKEND, DEVICE, ON_GPU
+
+torch.set_default_dtype(torch.float64)
+torch.set_default_device(DEVICE)
+
+
+def pytest_configure(config):
+    config.addinivalue_line("markers", "cpu_only: host-side test, skipped on a GPU")
+    if ON_GPU and not torch.cuda.is_available():
+        # a GPU run without a GPU must fail, not pass with every test on the CPU
+        pytest.exit(
+            f"MADSPACE_TEST_DEVICE={BACKEND}, but torch {torch.__version__} "
+            "has no GPU available",
+            returncode=1,
+        )
+
+
+def pytest_report_header(config):
+    name = torch.cuda.get_device_name(DEVICE) if ON_GPU else "cpu"
+    return f"madspace test device: {BACKEND} ({name})"
+
+
+def pytest_collection_modifyitems(config, items):
+    if not ON_GPU:
+        return
+    skip = pytest.mark.skip(reason=f"cpu_only, test device is {BACKEND}")
+    for item in items:
+        if "cpu_only" in item.keywords:
+            item.add_marker(skip)
diff --git a/madspace/tests/madspace_testing.py b/madspace/tests/madspace_testing.py
new file mode 100644
index 000000000..ccc5b2784
--- /dev/null
+++ b/madspace/tests/madspace_testing.py
@@ -0,0 +1,67 @@
+"""The device the tests run on, from MADSPACE_TEST_DEVICE: cpu (default), cuda or hip.
+
+torch creates its tensors on DEVICE by default (set in conftest.py), and madspace runs a
+function on the device of its inputs. torch calls the ROCm devices "cuda" too.
+"""
+
+import os
+
+import torch
+
+import madspace as ms
+
+BACKEND = os.environ.get("MADSPACE_TEST_DEVICE", "cpu").strip().lower() or "cpu"
+if BACKEND not in ("cpu", "cuda", "hip"):
+    raise ValueError(f"MADSPACE_TEST_DEVICE must be cpu, cuda or hip, not {BACKEND!r}")
+DEVICE = torch.device("cpu" if BACKEND == "cpu" else "cuda")
+ON_GPU = BACKEND != "cpu"
+
+
+def ms_device():
+    """The madspace device of the tests."""
+    return {"cpu": ms.cpu_device, "cuda": ms.cuda_device, "hip": ms.hip_device}[
+        BACKEND
+    ]()
+
+
+def default_context():
+    """The default madspace context of the device of the tests."""
+    return {
+        "cpu": ms.default_context,
+        "cuda": ms.default_cuda_context,
+        "hip": ms.default_hip_context,
+    }[BACKEND]()
+
+
+def new_context():
+    """A new madspace context on the device of the tests."""
+    return ms.Context(ms_device())
+
+
+def generator(seed):
+    """A seeded torch generator for rand, randn and randint below."""
+    return torch.Generator("cpu").manual_seed(seed)
+
+
+# The random numbers are drawn on the CPU and moved to the device of the tests, so that
+# every device gets the same inputs (and statistical tests the same outcome).
+def _on_device(draw, args, kwargs):
+    requires_grad = kwargs.pop("requires_grad", False)
+    return draw(*args, device="cpu", **kwargs).to(DEVICE).requires_grad_(requires_grad)
+
+
+def rand(*args, **kwargs):
+    return _on_device(torch.rand, args, kwargs)
+
+
+def randn(*args, **kwargs):
+    return _on_device(torch.randn, args, kwargs)
+
+
+def randint(*args, **kwargs):
+    return _on_device(torch.randint, args, kwargs)
+
+
+def to_numpy(tensor):
+    """A tensor as a numpy array on the host, e.g. for scipy."""
+    return tensor.detach().cpu().numpy()
diff --git a/madspace/tests/phasespace_helpers.py b/madspace/tests/phasespace_helpers.py
index ecceefe37..3f0a2732e 100644
--- a/madspace/tests/phasespace_helpers.py
+++ b/madspace/tests/phasespace_helpers.py
@@ -4,32 +4,31 @@
 __init__.py), so the test modules import it by name.
 """
 
-import numpy as np
+import torch
+from madspace_testing import generator, rand
 
 
 def sample(mapping, seed, n, conditions=()):
     """Map n uniform points through a PhaseSpaceMapping.
 
-    Returns the momenta and the weight (Jacobian times cuts) as arrays, as the
+    Returns the momenta and the weight (Jacobian times cuts) as tensors, as the
     mapping gives them: non-finite values are left in, see finite_weight.
     conditions are passed on unchanged, e.g. a permutation index per point.
     """
-    rng = np.random.default_rng(seed)
-    p_ext, _x1, _x2, det = mapping.map_forward(
-        [rng.random((n, mapping.random_dim()))], list(conditions)
-    )
-    return np.asarray(p_ext), np.asarray(det)
+    r = rand((n, mapping.random_dim()), generator=generator(seed))
+    p_ext, _x1, _x2, det = mapping.map_forward([r], list(conditions))
+    return p_ext, det
 
 
 def finite_weight(p, det):
     """The weight with every point whose weight or momenta are not finite set
     to 0, i.e. what an integral over the batch should sum."""
-    finite = np.isfinite(det) & np.all(np.isfinite(p), axis=(1, 2))
-    return np.where(finite, det, 0.0)
+    finite = torch.isfinite(det) & torch.isfinite(p).all(dim=2).all(dim=1)
+    return torch.where(finite, det, 0.0)
 
 
 def invariant_mass(p, *indices):
     """Invariant mass of the sum of the momenta at the given positions."""
     total = sum(p[:, i, :] for i in indices)
-    m2 = total[:, 0] ** 2 - np.sum(total[:, 1:] ** 2, axis=1)
-    return np.sqrt(np.maximum(m2, 0.0))
+    m2 = total[:, 0] ** 2 - torch.sum(total[:, 1:] ** 2, dim=1)
+    return torch.sqrt(torch.clamp(m2, min=0.0))
diff --git a/madspace/tests/test_2to3_scattering.py b/madspace/tests/test_2to3_scattering.py
index ed2244245..b72e39d98 100644
--- a/madspace/tests/test_2to3_scattering.py
+++ b/madspace/tests/test_2to3_scattering.py
@@ -1,7 +1,9 @@
+import math
 from collections import namedtuple
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand, randint, randn
 from pytest import approx
 
 import madspace as ms
@@ -13,27 +15,38 @@
 
 def inv_mass_sq(p):
     # p: (..., 4) with [E, px, py, pz]
-    return p[..., 0] ** 2 - np.sum(p[..., 1:] ** 2, axis=-1)
+    return p[..., 0] ** 2 - torch.sum(p[..., 1:] ** 2, dim=-1)
 
 
 def mass(p):
-    return np.sqrt(np.maximum(inv_mass_sq(p), 0.0))
+    return torch.sqrt(torch.clamp(inv_mass_sq(p), min=0.0))
 
 
-def boost(k: np.ndarray, p_boost: np.ndarray, inverse: bool = False) -> np.ndarray:
+def uniform(g, a, b, n):
+    return a + (b - a) * rand(n, generator=g)
+
+
+def assert_approx(actual, expected, rel=1e-6, abs=1e-12, msg=None):
+    """pytest.approx-like elementwise comparison of tensors."""
+    torch.testing.assert_close(actual, expected, rtol=rel, atol=abs, msg=msg)
+
+
+def boost(
+    k: torch.Tensor, p_boost: torch.Tensor, inverse: bool = False
+) -> torch.Tensor:
     # Change sign if inverse boost is performed
     sign = -1.0 if inverse else 1.0
     # Perform the boost
     EPS = 1e-16
-    rsq = np.clip(mass(p_boost), min=EPS)
-    kp = np.sum(k[..., 1:] * p_boost[..., 1:], axis=-1)
+    rsq = torch.clamp(mass(p_boost), min=EPS)
+    kp = torch.sum(k[..., 1:] * p_boost[..., 1:], dim=-1)
     k0 = (k[..., 0] * p_boost[..., 0] + sign * kp) / rsq
     c1 = (k[..., 0] + k0) / (rsq + p_boost[..., 0])
     k1 = k[..., 1] + sign * c1 * p_boost[..., 1]
     k2 = k[..., 2] + sign * c1 * p_boost[..., 2]
     k3 = k[..., 3] + sign * c1 * p_boost[..., 3]
 
-    return np.stack((k0, k1, k2, k3), axis=-1)
+    return torch.stack((k0, k1, k2, k3), dim=-1)
 
 
 # ----------------------------
@@ -61,15 +74,19 @@ def boost(k: np.ndarray, p_boost: np.ndarray, inverse: bool = False) -> np.ndarr
     ],
 )
 
-ZEROS = np.zeros(N)
-PA = np.full((N, 4), [np.sqrt(100**2 + 200**2 + 300**2), 100.0, 200.0, 300.0])
-PB = np.full((N, 4), [np.sqrt(100**2 + 200**2 + 400**2), 100.0, 200.0, -400.0])
+ZEROS = torch.zeros(N)
+PA = torch.tensor(
+    [math.sqrt(100**2 + 200**2 + 300**2), 100.0, 200.0, 300.0]
+).repeat(N, 1)
+PB = torch.tensor(
+    [math.sqrt(100**2 + 200**2 + 400**2), 100.0, 200.0, -400.0]
+).repeat(N, 1)
 P0 = PA + PB
 M0 = mass(P0)
-M1 = np.full(N, 10.0)
-M2 = np.full(N, 20.0)
-M3 = np.full(N, 30.0)
-M12 = np.full(N, 200.0)
+M1 = torch.full((N,), 10.0)
+M2 = torch.full((N,), 20.0)
+M3 = torch.full((N,), 30.0)
+M12 = torch.full((N,), 200.0)
 
 
 @pytest.fixture(
@@ -92,21 +109,21 @@ def fixed_input_points(rng, request):
     m0, m12, m1, m2, m3, pA, pB, p0 = request.param
 
     map_22 = ms.TwoToTwoParticleScattering(com=False)
-    r1 = rng.random(N)
-    r2 = rng.random(N)
+    r1 = rand(N, generator=rng)
+    r2 = rand(N, generator=rng)
     p12, p3, det_22 = map_22.map_forward([r1, r2, m12, m3], [pA, pB])
 
     # Discrete two-solution choice (int 0/1) + continuous randoms for the 2->3
-    r_choice = rng.integers(0, 2, size=N).astype(np.int32)
-    r_s23 = rng.random(N)
-    r_t1 = rng.random(N)
+    r_choice = randint(0, 2, (N,), generator=rng, dtype=torch.int32)
+    r_s23 = rand(N, generator=rng)
+    r_t1 = rand(N, generator=rng)
 
     return InputPoint(r_choice, r_s23, r_t1, m0, m12, m1, m2, m3, pA, pB, p0, p3, p12)
 
 
 @pytest.fixture
 def rng():
-    return np.random.default_rng(1234)
+    return generator(1234)
 
 
 @pytest.fixture(params=[False, True], ids=["base 2->2, LAB", "base 2->2, COM"])
@@ -119,40 +136,40 @@ def input_points(rng, request):
 
     # start with compatible incoming momenta
     # COM: pa=(E,0,0,+p), pb=(E,0,0,-p)
-    pz = rng.uniform(1000.0, 4000.0, N)
-    ma = rng.uniform(50.0, 300.0, N)
-    mb = rng.uniform(200.0, 500.0, N)
-    e1 = np.sqrt(pz**2 + ma**2)
-    e2 = np.sqrt(pz**2 + mb**2)
-    pa = np.stack([e1, np.zeros(N), np.zeros(N), +pz], axis=1)
-    pb = np.stack([e2, np.zeros(N), np.zeros(N), -pz], axis=1)
+    pz = uniform(rng, 1000.0, 4000.0, N)
+    ma = uniform(rng, 50.0, 300.0, N)
+    mb = uniform(rng, 200.0, 500.0, N)
+    e1 = torch.sqrt(pz**2 + ma**2)
+    e2 = torch.sqrt(pz**2 + mb**2)
+    pa = torch.stack([e1, torch.zeros(N), torch.zeros(N), +pz], dim=1)
+    pb = torch.stack([e2, torch.zeros(N), torch.zeros(N), -pz], dim=1)
 
     if not com:
         # boost to a some random frame
-        pa = np.array([[0, 0, 4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-        ea = np.sqrt(np.sum(pa**2, axis=1) + ma**2)
-        pa = np.concatenate([ea[:, None], pa], axis=1)
+        pa = torch.tensor([[0, 0, 4000.0]]) + 500.0 * randn((N, 3), generator=rng)
+        ea = torch.sqrt(torch.sum(pa**2, dim=1) + ma**2)
+        pa = torch.cat([ea[:, None], pa], dim=1)
 
-        pb = np.array([[0, 0, -4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-        eb = np.sqrt(np.sum(pb**2, axis=1) + mb**2)
-        pb = np.concatenate([eb[:, None], pb], axis=1)
+        pb = torch.tensor([[0, 0, -4000.0]]) + 500.0 * randn((N, 3), generator=rng)
+        eb = torch.sqrt(torch.sum(pb**2, dim=1) + mb**2)
+        pb = torch.cat([eb[:, None], pb], dim=1)
 
     # Build a valid 2->2 event to extract p3 as a spectator
     # Choose arbitrary outgoing masses for (Q, p3): make p3 light, Q heavyish
-    m3 = rng.uniform(1.0, 10.0, N)
-    m1 = rng.uniform(1.0, 40.0, N)
-    m2 = rng.uniform(1.0, 40.0, N)
-    m12 = rng.uniform(160, 400, N)
+    m3 = uniform(rng, 1.0, 10.0, N)
+    m1 = uniform(rng, 1.0, 40.0, N)
+    m2 = uniform(rng, 1.0, 40.0, N)
+    m12 = uniform(rng, 160, 400, N)
 
     map_22 = ms.TwoToTwoParticleScattering(com=com)
-    r1 = rng.random(N)
-    r2 = rng.random(N)
+    r1 = rand(N, generator=rng)
+    r2 = rand(N, generator=rng)
     p12, p3, det_22 = map_22.map_forward([r1, r2, m12, m3], [pa, pb])
 
     # Randoms for the 2->3 mapper
-    r_choice = rng.integers(0, 2, size=N).astype(np.int32)
-    r_s23 = rng.random(N)
-    r_t1 = rng.random(N)
+    r_choice = randint(0, 2, (N,), generator=rng, dtype=torch.int32)
+    r_s23 = rand(N, generator=rng)
+    r_t1 = rand(N, generator=rng)
 
     p0 = pa + pb
     m0 = mass(p0)
@@ -181,8 +198,8 @@ def test_momentum_conservation(input_points):
     p1, p2, det = mapping.map_forward(inputs, conditions)
     p_sum = p1 + p2 + input_points.p3
 
-    assert p_sum == approx(input_points.p0)
-    assert p1 + p2 == approx(input_points.p12)
+    assert_approx(p_sum, input_points.p0)
+    assert_approx(p1 + p2, input_points.p12)
 
 
 @pytest.mark.parametrize("arcsine", [False, True], ids=["flat s23", "arcsine s23"])
@@ -201,16 +218,14 @@ def test_inverse(input_points, arcsine):
     p1, p2, det = mapping.map_forward(inputs, conditions)
     *inv_inputs, inv_det = mapping.map_inverse([p1, p2], conditions)
 
-    assert inv_det == approx(1 / det, rel=1e-5)
+    assert_approx(inv_det, 1 / det, rel=1e-5)
 
     for i, (inp, inv_inp) in enumerate(zip(inputs, inv_inputs)):
         if i == 0:
             # discrete two-solution choice: recovered exactly as an int
-            assert (
-                np.asarray(inv_inp).astype(np.int64) == np.asarray(inp).astype(np.int64)
-            ).all()
+            assert torch.equal(inv_inp.long(), inp.long())
             continue
-        assert inp == approx(inv_inp), f"mismatch in input index {i}"
+        assert_approx(inp, inv_inp, msg=f"mismatch in input index {i}")
 
 
 def test_on_shell_masses(input_points):
@@ -228,16 +243,16 @@ def test_on_shell_masses(input_points):
     p1, p2, det = mapping.map_forward(inputs, conditions)
 
     # Outgoing masses must match m1, m2; spectator stays whatever it was.
-    assert mass(p1 + p2) == approx(input_points.m12)
-    assert mass(p1) == approx(input_points.m1)
-    assert mass(p2) == approx(input_points.m2)
+    assert_approx(mass(p1 + p2), input_points.m12)
+    assert_approx(mass(p1), input_points.m1)
+    assert_approx(mass(p2), input_points.m2)
 
 
 def test_phase_space_compare(rng, input_points):
     mapping23 = ms.TwoToThreeParticleScattering()
     mapping22 = ms.TwoToTwoParticleScattering(com=False)
-    r1 = rng.random(N)
-    r2 = rng.random(N)
+    r1 = rand(N, generator=rng)
+    r2 = rand(N, generator=rng)
 
     inputs = [
         input_points.r_choice,
@@ -263,8 +278,10 @@ def test_phase_space_compare(rng, input_points):
     det23 = det23 * 2.0
 
     # Outgoing masses must match m1, m2; spectator stays whatever it was.
-    std_error_23 = np.std(det23) / np.sqrt(N)
-    assert np.mean(det23) == approx(np.mean(det22), abs=3 * std_error_23, rel=1e-6)
+    std_error_23 = det23.std(correction=0).item() / math.sqrt(N)
+    assert det23.mean().item() == approx(
+        det22.mean().item(), abs=3 * std_error_23, rel=1e-6
+    )
 
 
 def test_phase_space_volume(fixed_input_points):
@@ -288,7 +305,10 @@ def test_phase_space_volume(fixed_input_points):
     m1_2 = fixed_input_points.m1**2
     m2_2 = fixed_input_points.m2**2
     phase_space_vol = (
-        np.pi / (2 * s) * np.sqrt((s - m1_2 - m2_2) ** 2 - 4 * m1_2 * m2_2)
+        math.pi / (2 * s) * torch.sqrt((s - m1_2 - m2_2) ** 2 - 4 * m1_2 * m2_2)
+    )
+    std_error = det.std(correction=0).item() / math.sqrt(N)
+    # fixed masses: the volume is the same for every point
+    assert det.mean().item() == approx(
+        phase_space_vol[0].item(), abs=3 * std_error, rel=1e-6
     )
-    std_error = np.std(det) / np.sqrt(N)
-    assert np.mean(det) == approx(phase_space_vol, abs=3 * std_error, rel=1e-6)
diff --git a/madspace/tests/test_buffer_unweighter.py b/madspace/tests/test_buffer_unweighter.py
index 5ae349d67..514c7b87e 100644
--- a/madspace/tests/test_buffer_unweighter.py
+++ b/madspace/tests/test_buffer_unweighter.py
@@ -20,7 +20,8 @@
 binomial, and checked to 5 sigma.
 """
 
-import numpy as np
+import torch
+from madspace_testing import default_context, generator, rand, randn
 
 import madspace as ms
 
@@ -37,19 +38,19 @@ def run(weights, prob, quantile=QUANTILE):
         ]
     )
     runtime = ms.FunctionRuntime(
-        ms.BufferUnweighter(types, quantile).function(), ms.default_context()
+        ms.BufferUnweighter(types, quantile).function(), default_context()
     )
-    index = np.arange(len(weights), dtype=np.float64)
+    index = torch.arange(len(weights), dtype=torch.float64)
     uw, prob_out, index_out = (
-        ms.Tensor.numpy(t) for t in runtime.call([weights, prob, index])
+        torch.from_dlpack(t) for t in runtime.call([weights, prob, index])
     )
-    return uw, prob_out, index_out.astype(np.int64)
+    return uw, prob_out, index_out.long()
 
 
 def expected_max_weight(weights, quantile=QUANTILE):
     """What op_quantile picks: element floor(q N) of the sorted |w|."""
     position = min(int(quantile * len(weights)), len(weights) - 1)
-    return np.sort(np.abs(weights))[position]
+    return torch.sort(weights.abs()).values[position]
 
 
 def check(weights, prob):
@@ -61,43 +62,43 @@ def check(weights, prob):
     assert w_max > 0
 
     # the sign is kept, the weight is max(|w|, w_max), and w_max is reached
-    np.testing.assert_array_equal(np.sign(uw), np.sign(w))
-    np.testing.assert_array_equal(np.abs(uw), np.maximum(np.abs(w), w_max))
-    assert np.abs(uw).min() == w_max
+    assert torch.equal(torch.sign(uw), torch.sign(w))
+    assert torch.equal(uw.abs(), torch.maximum(w.abs(), w_max))
+    assert uw.abs().min() == w_max
 
     # the density of the kept events: acceptance probability times N / N_kept
-    np.testing.assert_allclose(
-        prob_out, w / uw * q * (len(weights) / len(index)), rtol=1e-14
+    torch.testing.assert_close(
+        prob_out, w / uw * q * (len(weights) / len(index)), rtol=1e-14, atol=0
     )
 
     # the kept count is a sum of Bernoulli draws with p = min(1, |w| / w_max)
-    p = np.minimum(1.0, np.abs(weights) / w_max)
-    assert abs(len(index) - p.sum()) < 5 * np.sqrt(np.sum(p * (1 - p))) + 1
+    p = torch.clamp(weights.abs() / w_max, max=1.0)
+    assert abs(len(index) - p.sum()) < 5 * torch.sqrt(torch.sum(p * (1 - p))) + 1
 
     # sum_kept uw / N estimates the mean weight (each event contributes
     # sign(w) max(|w|, w_max) with probability p, i.e. w on average)
-    spread = np.sqrt(np.sum((1 - p) * np.abs(weights) * w_max)) / len(weights)
+    spread = torch.sqrt(torch.sum((1 - p) * weights.abs() * w_max)) / len(weights)
     assert abs(uw.sum() / len(weights) - weights.mean()) < 5 * spread
     return len(index) / len(weights)
 
 
-def batch(rng, scale=1.0):
-    return scale * rng.lognormal(0.0, 1.0, N)
+def batch(g, scale=1.0):
+    return scale * torch.exp(randn(N, generator=g))
 
 
 def test_positive_batch():
-    rng = np.random.default_rng(1)
-    weights = batch(rng)
-    kept = check(weights, rng.uniform(0.5, 1.5, N))
+    g = generator(1)
+    weights = batch(g)
+    kept = check(weights, 0.5 + rand(N, generator=g))
     assert kept < 0.5
 
 
 def test_negative_batch_is_unweighted():
     """A channel where the interference is negative throughout: the signed
     quantile was negative, so every event was kept at its full weight."""
-    rng = np.random.default_rng(2)
-    weights = -batch(rng)
-    kept = check(weights, rng.uniform(0.5, 1.5, N))
+    g = generator(2)
+    weights = -batch(g)
+    kept = check(weights, 0.5 + rand(N, generator=g))
     assert kept < 0.5
 
 
@@ -105,8 +106,8 @@ def test_mixed_batch_is_unweighted():
     """85% large negative weights and a short positive tail: the signed 95%
     quantile came from the positive tail, about 30 times below the |w| one,
     and nearly every negative event was kept as an over-weight."""
-    rng = np.random.default_rng(3)
-    negative = rng.uniform(size=N) < 0.85
-    weights = np.where(negative, -batch(rng, 10.0), batch(rng, 1.0))
-    kept = check(weights, rng.uniform(0.5, 1.5, N))
+    g = generator(3)
+    negative = rand(N, generator=g) < 0.85
+    weights = torch.where(negative, -batch(g, 10.0), batch(g, 1.0))
+    kept = check(weights, 0.5 + rand(N, generator=g))
     assert kept < 0.5
diff --git a/madspace/tests/test_color_ordered_pair_cuts.py b/madspace/tests/test_color_ordered_pair_cuts.py
index d85c15520..f2d58a3de 100644
--- a/madspace/tests/test_color_ordered_pair_cuts.py
+++ b/madspace/tests/test_color_ordered_pair_cuts.py
@@ -15,8 +15,9 @@
 
 import math
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand, randint
 
 import madspace as ms
 
@@ -65,9 +66,9 @@ def matrix(pairs, value):
 
 
 def random_inputs(random_dim, discrete_dim, n, seed):
-    rng = np.random.default_rng(seed)
-    r = rng.random((n, random_dim))
-    d = rng.integers(0, 2, size=(n, discrete_dim)).astype(np.int32)
+    g = generator(seed)
+    r = rand((n, random_dim), generator=g)
+    d = randint(0, 2, (n, discrete_dim), generator=g, dtype=torch.int32)
     return r, d
 
 
@@ -83,17 +84,17 @@ def psm_momenta(order, cut_items, r, d):
     assert mapping.random_dim() == r.shape[1]
     assert mapping.discrete_dim() == d.shape[1]
     inputs = [r, d] if d.shape[1] else [r]
-    return np.asarray(mapping.map_forward(inputs)[0])
+    return mapping.map_forward(inputs)[0]
 
 
 def co_momenta(order, pt_min, m_inv_min, dr_min, r, d):
     mapping = ms.ColorOrderedMapping(order, 0.8, 0.8, pt_min, m_inv_min, dr_min)
     n = r.shape[0]
-    inputs = [r[:, i].copy() for i in range(r.shape[1])]
-    inputs += [d[:, j].copy() for j in range(d.shape[1])]
-    conditions = [np.full(n, CM_ENERGY)] + [np.zeros(n) for _ in range(N_OUT)]
+    inputs = [r[:, i].contiguous() for i in range(r.shape[1])]
+    inputs += [d[:, j].contiguous() for j in range(d.shape[1])]
+    conditions = [torch.full((n,), CM_ENERGY)] + [torch.zeros(n) for _ in range(N_OUT)]
     out = mapping.map_forward(inputs, conditions)
-    return np.stack([np.asarray(p) for p in out[:-1]], axis=1)
+    return torch.stack(out[:-1], dim=1)
 
 
 def same_points(p, q):
@@ -101,9 +102,9 @@ def same_points(p, q):
 
     The bound narrows ranges, so some points land where nothing is left to
     sample and come out as NaN (weight zero); those count as the same point."""
-    nan_p = ~np.all(np.isfinite(p), axis=(1, 2))
-    nan_q = ~np.all(np.isfinite(q), axis=(1, 2))
-    close = np.all(np.abs(np.nan_to_num(p - q)) <= ATOL, axis=(1, 2))
+    nan_p = ~torch.isfinite(p).all(dim=2).all(dim=1)
+    nan_q = ~torch.isfinite(q).all(dim=2).all(dim=1)
+    close = (torch.nan_to_num(p - q).abs() <= ATOL).all(dim=2).all(dim=1)
     return (nan_p & nan_q) | (~nan_p & ~nan_q & close)
 
 
@@ -146,13 +147,13 @@ def test_pair_cut_reaches_color_ordered_mapping(
 
     # Near degenerate configurations the extra operations in
     # PhaseSpaceMapping can be amplified past ATOL, so allow a few points.
-    assert np.mean(same_points(p_psm, p_with)) > 0.99
+    assert same_points(p_psm, p_with).double().mean() > 0.99
 
     # Without the comparison below the test would pass for a cut that does
     # not affect this topology at all, and prove nothing. A floor only moves
     # the points where it lies above the kinematic one: all of them for the
     # pair masses here, about half for the delta-R cut.
-    assert np.mean(same_points(p_with, p_without)) < 0.9
+    assert same_points(p_with, p_without).double().mean() < 0.9
 
 
 def sample(order, cuts, n, seed):
@@ -167,14 +168,14 @@ def sample(order, cuts, n, seed):
     r, d = random_inputs(mapping.random_dim(), mapping.discrete_dim(), n, seed)
     inputs = [r, d] if mapping.discrete_dim() else [r]
     p_ext, _x1, _x2, det = mapping.map_forward(inputs)
-    det = np.nan_to_num(np.asarray(det)) * 2.0 ** mapping.discrete_dim()
-    return np.asarray(p_ext), det
+    det = torch.nan_to_num(det) * 2.0 ** mapping.discrete_dim()
+    return p_ext, det
 
 
 def invariant_mass(p, a, b):
     total = p[:, a + 2] + p[:, b + 2]
-    m2 = total[:, 0] ** 2 - np.sum(total[:, 1:] ** 2, axis=1)
-    return np.sqrt(np.maximum(m2, 0.0))
+    m2 = total[:, 0] ** 2 - torch.sum(total[:, 1:] ** 2, dim=1)
+    return torch.sqrt(torch.clamp(m2, min=0.0))
 
 
 # The central 2->2 split, whose set masses are the only place these cuts
@@ -188,17 +189,17 @@ def test_pair_cut_on_odd_child_raises_efficiency_not_volume(pair):
     a, b = pair
 
     p_free, det_free = sample(SPLIT_2_2, None, n, SEED)
-    external = np.where(invariant_mass(p_free, a, b) >= PAIR_MASS_CUT, det_free, 0.0)
+    external = torch.where(invariant_mass(p_free, a, b) >= PAIR_MASS_CUT, det_free, 0.0)
 
     p_cut, piped = sample(SPLIT_2_2, ms.Cuts([pair_mass(a, b)]), n, SEED + 1)
-    assert np.all(invariant_mass(p_cut, a, b)[piped != 0.0] >= PAIR_MASS_CUT)
-    assert np.all(piped >= 0.0)
+    assert torch.all(invariant_mass(p_cut, a, b)[piped != 0.0] >= PAIR_MASS_CUT)
+    assert torch.all(piped >= 0.0)
 
     # with the bound dropped, both fractions are the same (0.84 / 0.17)
-    eff_external = np.mean(external > 0.0)
-    eff_piped = np.mean(piped > 0.0)
+    eff_external = (external > 0.0).double().mean()
+    eff_piped = (piped > 0.0).double().mean()
     assert eff_piped > eff_external + 0.03
 
-    error = math.sqrt((external.var() + piped.var()) / n)
+    error = math.sqrt((external.var(correction=0) + piped.var(correction=0)) / n)
     assert error > 0.0
-    assert abs(piped.mean() - external.mean()) < 5.0 * error
+    assert abs(piped.mean() - external.mean()).item() < 5.0 * error
diff --git a/madspace/tests/test_color_ordered_precision.py b/madspace/tests/test_color_ordered_precision.py
index bab5d19ab..63e4b2dc5 100644
--- a/madspace/tests/test_color_ordered_precision.py
+++ b/madspace/tests/test_color_ordered_precision.py
@@ -10,8 +10,9 @@
   inverse reads a massless particle's mass off its momentum.
 """
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand, randint
 
 import madspace as ms
 
@@ -27,14 +28,14 @@ def run(color_order, cm_energy, seed, cuts=None):
     cuts = cuts or {}
     n_out = len(color_order) - 2
     mapping = ms.ColorOrderedMapping(color_order, **cuts)
-    rng = np.random.default_rng(seed)
-    r = rng.random((N, mapping.random_dim()))
-    d = rng.integers(0, 2, size=(N, mapping.discrete_dim())).astype(np.int32)
-    inputs = [r[:, i].copy() for i in range(r.shape[1])]
-    inputs += [d[:, j].copy() for j in range(d.shape[1])]
-    conditions = [np.full(N, cm_energy)] + [np.zeros(N)] * n_out
+    g = generator(seed)
+    r = rand((N, mapping.random_dim()), generator=g)
+    d = randint(0, 2, (N, mapping.discrete_dim()), generator=g, dtype=torch.int32)
+    inputs = [r[:, i].contiguous() for i in range(r.shape[1])]
+    inputs += [d[:, j].contiguous() for j in range(d.shape[1])]
+    conditions = [torch.full((N,), cm_energy)] + [torch.zeros(N)] * n_out
     *momenta, det = mapping.map_forward(inputs, conditions)
-    return mapping, inputs, conditions, momenta, r, np.asarray(det)
+    return mapping, inputs, conditions, momenta, r, det
 
 
 @pytest.mark.parametrize(
@@ -45,29 +46,30 @@ def run(color_order, cm_energy, seed, cuts=None):
 def test_massless_on_shell(color_order):
     cm_energy = 13000.0
     *_, momenta, _, det = run(color_order, cm_energy, 5)
-    p = np.stack([np.asarray(k) for k in momenta], axis=1)[:, 2:]
+    p = torch.stack(momenta, dim=1)[:, 2:]
     p = p[det > 0]
-    m2 = p[..., 0] ** 2 - np.sum(p[..., 1:] ** 2, axis=-1)
+    m2 = p[..., 0] ** 2 - torch.sum(p[..., 1:] ** 2, dim=-1)
     # sqrt(|m^2|) / E is ~1e-8 from rounding alone; it reached 4e-5 before
-    assert np.max(np.sqrt(np.abs(m2)) / p[..., 0]) < 1e-6
+    assert torch.max(torch.sqrt(torch.abs(m2)) / p[..., 0]) < 1e-6
     # momentum conservation to a few ulps of the total energy
-    p_sum = np.sum(p, axis=1)
+    p_sum = torch.sum(p, dim=1)
     p_sum[:, 0] -= cm_energy
-    assert np.max(np.abs(p_sum)) < 1e-14 * cm_energy
+    assert torch.max(torch.abs(p_sum)) < 1e-14 * cm_energy
 
 
 def test_round_trip_with_cuts():
     order = [0, 2, 3, 4, 5, 1]
     mapping, inputs, conditions, momenta, r, det = run(order, 1000.0, 11, CUTS)
     out = mapping.map_inverse(momenta, conditions)
-    r_inv = np.stack([np.asarray(o) for o in out[: mapping.random_dim()]], axis=1)
-    det_inv = np.asarray(out[-1])
+    r_inv = torch.stack(out[: mapping.random_dim()], dim=1)
+    det_inv = out[-1]
     ok = det > 0
-    rt = np.abs(det * det_inv - 1.0)[ok]
+    rt = torch.abs(det * det_inv - 1.0)[ok]
     # about 3% of the points failed before, carrying ~4% of the fiducial volume
-    assert np.mean(rt > 1e-2) < 1e-4
-    assert np.quantile(rt, 0.999) < 1e-6
-    assert np.mean(np.max(np.abs(r_inv - r), axis=1)[ok] > 1e-6) < 1e-4
+    assert (rt > 1e-2).double().mean() < 1e-4
+    assert torch.quantile(rt, 0.999) < 1e-6
+    r_err = torch.max(torch.abs(r_inv - r), dim=1).values[ok]
+    assert (r_err > 1e-6).double().mean() < 1e-4
 
 
 @pytest.mark.parametrize(
@@ -87,19 +89,19 @@ def test_lab_frame_on_shell(masses):
         mode=ms.PhaseSpaceMapping.color_ordered,
         color_order=[0, 2, 3, 4, 5, 1],
     )
-    rng = np.random.default_rng(3)
-    inputs = [rng.random((N, mapping.random_dim()))]
+    g = generator(3)
+    inputs = [rand((N, mapping.random_dim()), generator=g)]
     if mapping.discrete_dim():
         inputs.append(
-            rng.integers(0, 2, size=(N, mapping.discrete_dim())).astype(np.int32)
+            randint(0, 2, (N, mapping.discrete_dim()), generator=g, dtype=torch.int32)
         )
     p, x1, x2, det = mapping.map_forward(inputs)
     # zero-weight points at x1 or x2 ~ 1e-13 have no meaningful momenta
     p = p[det > 1e-20]
-    m2 = p[..., 0] ** 2 - np.sum(p[..., 1:] ** 2, axis=-1)
+    m2 = p[..., 0] ** 2 - torch.sum(p[..., 1:] ** 2, dim=-1)
     # |p^2 - m^2| / E^2: a few 1e-15 for massless momenta, ~1e-12 for the top
     # that absorbs the rounding of the others; it reached ~1e-9 before
-    off_shell = np.abs(m2 - np.square(masses)) / p[..., 0] ** 2
-    assert np.max(off_shell) < 1e-11
-    violation = np.abs(np.sum(p[:, 2:], axis=1) - np.sum(p[:, :2], axis=1))
-    assert np.max(violation) < 1e-14 * cm_energy
+    off_shell = torch.abs(m2 - torch.tensor(masses) ** 2) / p[..., 0] ** 2
+    assert torch.max(off_shell) < 1e-11
+    violation = torch.abs(torch.sum(p[:, 2:], dim=1) - torch.sum(p[:, :2], dim=1))
+    assert torch.max(violation) < 1e-14 * cm_energy
diff --git a/madspace/tests/test_cut_kernels.py b/madspace/tests/test_cut_kernels.py
index 24ad600bf..b605cd990 100644
--- a/madspace/tests/test_cut_kernels.py
+++ b/madspace/tests/test_cut_kernels.py
@@ -15,8 +15,9 @@
 import os
 import re
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import default_context
 
 import madspace as ms
 
@@ -25,7 +26,7 @@
 
 # rows: all inside, one outside below, one outside above, all outside, on the
 # edges (inside: the window is closed)
-OBS = np.array(
+OBS = torch.tensor(
     [
         [1.2, 1.5, 1.9, 1.1],
         [0.5, 1.5, 1.9, 1.1],
@@ -44,17 +45,18 @@ def run_cuts():
     obs = fb.input(0)
     fb.output(0, fb.cut_all(obs, LOW, HIGH))
     fb.output(1, fb.cut_any(obs, LOW, HIGH))
-    runtime = ms.FunctionRuntime(fb.function(), ms.default_context())
-    return [ms.Tensor.numpy(t) for t in runtime.call([OBS])]
+    runtime = ms.FunctionRuntime(fb.function(), default_context())
+    return [torch.from_dlpack(t) for t in runtime.call([OBS])]
 
 
 def test_cut_all_and_any_apply_one_window_to_every_observable():
     cut_all, cut_any = run_cuts()
     inside = (OBS >= LOW) & (OBS <= HIGH)
-    np.testing.assert_array_equal(cut_all, inside.all(axis=1).astype(float))
-    np.testing.assert_array_equal(cut_any, inside.any(axis=1).astype(float))
+    torch.testing.assert_close(cut_all, inside.all(dim=1).double(), rtol=0, atol=0)
+    torch.testing.assert_close(cut_any, inside.any(dim=1).double(), rtol=0, atol=0)
 
 
+@pytest.mark.cpu_only
 def test_kernel_inputs_have_the_declared_dimensions():
     yaml = pytest.importorskip("yaml")
     root = os.path.join(os.path.dirname(os.path.abspath(__file__)), os.pardir)
diff --git a/madspace/tests/test_cut_rapidity_windows.py b/madspace/tests/test_cut_rapidity_windows.py
index 6f2ec0dbe..46554d0f5 100644
--- a/madspace/tests/test_cut_rapidity_windows.py
+++ b/madspace/tests/test_cut_rapidity_windows.py
@@ -25,9 +25,9 @@
 
 import math
 
-import numpy as np
 import pytest
-from pytest import approx
+import torch
+from madspace_testing import generator, rand, randint
 
 import madspace as ms
 
@@ -93,46 +93,64 @@ def jet_cuts(pids, pt=PT_JET, eta=ETA_JET):
     return items
 
 
+def finite(p_ext):
+    return torch.isfinite(p_ext).all(dim=2).all(dim=1)
+
+
 def sample(mapping, seed=SEED, n=BATCH_SIZE):
-    rng = np.random.default_rng(seed)
-    r = rng.random((n, mapping.random_dim()))
+    r = rand((n, mapping.random_dim()), generator=generator(seed))
     p_ext, x1, x2, det = mapping.map_forward([r])
-    p_ext, x1, x2, det = (np.asarray(a) for a in (p_ext, x1, x2, det))
-    ok = np.isfinite(det) & np.all(np.isfinite(p_ext), axis=(1, 2))
-    return r, p_ext, x1, x2, np.where(ok, det, 0.0)
+    ok = torch.isfinite(det) & finite(p_ext)
+    return r, p_ext, x1, x2, torch.where(ok, det, 0.0)
 
 
 def passes(cuts, p_ext):
-    return np.asarray(cuts(p_ext)).reshape(-1) > 0.5
+    return cuts(p_ext).reshape(-1) > 0.5
+
+
+def fraction(mask):
+    """The fraction of True entries of a boolean tensor."""
+    return mask.double().mean().item()
+
+
+def mean_and_error(w):
+    """Mean of a batch and its standard error (numpy's std, without Bessel)."""
+    return w.mean().item(), w.std(correction=0).item() / math.sqrt(len(w))
 
 
 def pt(p):
-    return np.hypot(p[..., 1], p[..., 2])
+    return torch.hypot(p[..., 1], p[..., 2])
 
 
 def rapidity(p):
-    return 0.5 * np.log((p[..., 0] + p[..., 3]) / (p[..., 0] - p[..., 3]))
+    return 0.5 * torch.log((p[..., 0] + p[..., 3]) / (p[..., 0] - p[..., 3]))
 
 
 def eta(p):
-    p_mag = np.sqrt(np.sum(p[..., 1:] ** 2, axis=-1))
-    return 0.5 * np.log((p_mag + p[..., 3]) / (p_mag - p[..., 3]))
+    p_mag = torch.sqrt(torch.sum(p[..., 1:] ** 2, dim=-1))
+    return 0.5 * torch.log((p_mag + p[..., 3]) / (p_mag - p[..., 3]))
 
 
 def beam_rapidity(x1, x2):
-    return 0.5 * np.log(x1 / x2)
+    return 0.5 * torch.log(x1 / x2)
 
 
 def sqrt_s_hat(x1, x2):
-    return CM_ENERGY * np.sqrt(x1 * x2)
+    return CM_ENERGY * torch.sqrt(x1 * x2)
+
+
+def assert_close(actual, expected, rel=0.0, abs=0.0):
+    """Like pytest.approx for tensors: |actual - expected| <= abs + rel * |expected|."""
+    expected = torch.as_tensor(expected, dtype=actual.dtype).expand_as(actual)
+    torch.testing.assert_close(actual, expected, rtol=rel, atol=abs)
 
 
 def weighted_cut_integral(mapping, cuts, seed, n=400_000):
     """det / s_hat over the cut region: a propagator-shaped test integrand,
     so the low-mass region carries weight."""
     _, p_ext, x1, x2, det = sample(mapping, seed=seed, n=n)
-    weight = np.where(passes(cuts, p_ext), det / (CM_ENERGY**2 * x1 * x2), 0.0)
-    return weight.mean(), weight.std() / math.sqrt(n)
+    weight = torch.where(passes(cuts, p_ext), det / (CM_ENERGY**2 * x1 * x2), 0.0)
+    return mean_and_error(weight)
 
 
 def assert_same_integral(topology, cuts, n=400_000):
@@ -160,13 +178,15 @@ def test_dy_every_generated_point_passes():
     windows, so with the cuts handed to the mapping every point it generates
     with a weight passes them, while without them most fail."""
     c = dy_cuts()
-    _, p_free, _, _, det_free = sample(ms.PhaseSpaceMapping(s_channel_two_body(), CM_ENERGY))
-    assert np.mean(passes(c, p_free)[det_free > 0]) < 0.9
+    _, p_free, _, _, det_free = sample(
+        ms.PhaseSpaceMapping(s_channel_two_body(), CM_ENERGY)
+    )
+    assert fraction(passes(c, p_free)[det_free > 0]) < 0.9
     _, p_ext, _, _, det = sample(
         ms.PhaseSpaceMapping(s_channel_two_body(), CM_ENERGY, cuts=c)
     )
-    assert np.mean(det > 0) > 0.999
-    assert np.mean(passes(c, p_ext)[det > 0]) > 0.999
+    assert fraction(det > 0) > 0.999
+    assert fraction(passes(c, p_ext)[det > 0]) > 0.999
 
 
 def test_dy_windows_are_reached():
@@ -178,10 +198,10 @@ def test_dy_windows_are_reached():
     )
     leptons = p_ext[:, 2:4]
     y = beam_rapidity(x1, x2)
-    assert np.abs(y).max() <= ETA_LEPTON + 1e-9
-    assert np.abs(y).max() > 0.99 * ETA_LEPTON
-    assert np.abs(eta(leptons)).max() <= ETA_LEPTON + 1e-6
-    assert np.abs(eta(leptons)).max() > 0.999 * ETA_LEPTON
+    assert torch.abs(y).max() <= ETA_LEPTON + 1e-9
+    assert torch.abs(y).max() > 0.99 * ETA_LEPTON
+    assert torch.abs(eta(leptons)).max() <= ETA_LEPTON + 1e-6
+    assert torch.abs(eta(leptons)).max() > 0.999 * ETA_LEPTON
     assert pt(leptons).min() >= PT_LEPTON * (1 - 1e-9)
     assert pt(leptons).min() < 1.01 * PT_LEPTON
 
@@ -195,10 +215,9 @@ def test_dy_round_trip():
     r, p_ext, x1, x2, det = sample(mapping, n=10_000)
     keep = det > 0
     r_back, det_back = mapping.map_inverse([p_ext, x1, x2], [])
-    r_back, det_back = np.asarray(r_back), np.asarray(det_back)
     # the cut mask is part of det but not of det_back
-    assert r_back[keep] == approx(r[keep], abs=1e-7)
-    assert (det * det_back)[keep] == approx(1.0, rel=1e-7)
+    assert_close(r_back[keep], r[keep], abs=1e-7)
+    assert_close((det * det_back)[keep], 1.0, rel=1e-7)
 
 
 # --------------------------------------------------------------------------
@@ -215,12 +234,12 @@ def test_massive_partner_without_cuts(pt_jet):
     c = ms.Cuts(jet_cuts(ZG_PIDS, pt=pt_jet))
     topology = s_channel_two_body((M_Z, 0.0))
     _, p_free, _, _, _ = sample(ms.PhaseSpaceMapping(topology, CM_ENERGY))
-    assert np.mean(passes(c, p_free)) < 0.9
+    assert fraction(passes(c, p_free)) < 0.9
     _, p_ext, x1, x2, det = sample(ms.PhaseSpaceMapping(topology, CM_ENERGY, cuts=c))
-    assert np.mean(passes(c, p_ext)) > 0.99
-    assert np.mean(passes(c, p_ext)[det > 0]) > 0.999
+    assert fraction(passes(c, p_ext)) > 0.99
+    assert fraction(passes(c, p_ext)[det > 0]) > 0.999
     # Y is not restricted
-    assert np.abs(beam_rapidity(x1, x2)[det > 0]).max() > ETA_JET + 0.5
+    assert torch.abs(beam_rapidity(x1, x2)[det > 0]).max() > ETA_JET + 0.5
     assert_same_integral(topology, c)
 
 
@@ -272,12 +291,12 @@ def test_composite_rapidity_windows():
     topology = s_channel_four_body()
     _, p_ext, x1, x2, det = sample(ms.PhaseSpaceMapping(topology, CM_ENERGY, cuts=c))
     # every sampled point: the weight already contains the cut mask
-    keep = np.all(np.isfinite(p_ext), axis=(1, 2))
+    keep = finite(p_ext)
     eta_max = max(ETA_LEPTON, ETA_JET)
-    assert np.abs(beam_rapidity(x1, x2)[keep]).max() <= eta_max + 1e-9
+    assert torch.abs(beam_rapidity(x1, x2)[keep]).max() <= eta_max + 1e-9
     for pair in ((2, 3), (4, 5)):
         system = p_ext[keep, pair[0]] + p_ext[keep, pair[1]]
-        assert np.abs(rapidity(system)).max() <= eta_max + 1e-6
+        assert torch.abs(rapidity(system)).max() <= eta_max + 1e-6
     assert_same_integral(topology, c)
 
 
@@ -286,13 +305,14 @@ def test_composite_round_trip():
     r, p_ext, x1, x2, det = sample(mapping, n=10_000)
     keep = det > 0
     r_back, det_back = mapping.map_inverse([p_ext, x1, x2], [])
-    r_back, det_back = np.asarray(r_back), np.asarray(det_back)
-    # a handful of points lose ~1e-5 in the inner decays with or without the
-    # windows (momentum reconstruction), hence quantiles
-    r_diff = np.abs(r_back[keep] - r[keep])
-    det_diff = np.abs((det * det_back)[keep] - 1.0)
-    assert np.quantile(r_diff, 0.999) < 1e-7 and r_diff.max() < 1e-4
-    assert np.quantile(det_diff, 0.99) < 1e-7 and np.quantile(det_diff, 0.999) < 1e-5
+    # a handful of points lose ~1e-5 to 1e-4 (platform dependent) in the inner decays
+    # with or without the windows (momentum reconstruction), hence quantiles
+    r_diff = torch.abs(r_back[keep] - r[keep])
+    det_diff = torch.abs((det * det_back)[keep] - 1.0)
+    assert torch.quantile(r_diff, 0.999) < 1e-7 and r_diff.max() < 1e-3
+    assert (
+        torch.quantile(det_diff, 0.99) < 1e-7 and torch.quantile(det_diff, 0.999) < 1e-5
+    )
 
 
 # --------------------------------------------------------------------------
@@ -308,9 +328,9 @@ def test_no_rapidity_window_with_an_uncut_particle():
         ms.PhaseSpaceMapping(s_channel_four_body(), CM_ENERGY, cuts=c)
     )
     keep = det > 0
-    assert np.abs(beam_rapidity(x1, x2)[keep]).max() > ETA_LEPTON + 1.0
+    assert torch.abs(beam_rapidity(x1, x2)[keep]).max() > ETA_LEPTON + 1.0
     gluons = p_ext[keep, 4] + p_ext[keep, 5]
-    assert np.abs(rapidity(gluons)).max() > ETA_LEPTON + 1.0
+    assert torch.abs(rapidity(gluons)).max() > ETA_LEPTON + 1.0
 
 
 def test_no_cuts_no_change():
@@ -319,7 +339,7 @@ def test_no_cuts_no_change():
     a = sample(ms.PhaseSpaceMapping(topology, CM_ENERGY))
     b = sample(ms.PhaseSpaceMapping(topology, CM_ENERGY, cuts=ms.Cuts(4)))
     for u, v in zip(a, b):
-        assert np.array_equal(u, v)
+        assert torch.equal(u, v)
 
 
 def test_two_body_decay_window_alone():
@@ -329,21 +349,20 @@ def test_two_body_decay_window_alone():
     pt_min = 30.0
     decay = ms.TwoBodyDecay(True, pt_min=pt_min)
     n = 10_000
-    rng = np.random.default_rng(SEED)
-    r_phi, r_cos = rng.random(n), rng.random(n)
-    m0 = np.full(n, 100.0)
-    zeros = np.zeros(n)
-    ones = np.ones(n)
+    g = generator(SEED)
+    r_phi, r_cos = rand(n, generator=g), rand(n, generator=g)
+    m0 = torch.full((n,), 100.0)
+    zeros = torch.zeros(n)
+    ones = torch.ones(n)
     p1, p2, det = decay.map_forward([r_phi, r_cos, m0, zeros, zeros], [ones, ones])
-    p1, det = np.asarray(p1), np.asarray(det)
     plain = ms.TwoBodyDecay(True)
     _, _, det_plain = plain.map_forward([r_phi, r_cos, m0, zeros, zeros], [])
     kept = math.sqrt(1 - (pt_min / 50.0) ** 2)
-    assert det == approx(np.asarray(det_plain) * kept, rel=1e-12)
+    assert_close(det, det_plain * kept, rel=1e-12)
     assert pt(p1).min() >= pt_min * (1 - 1e-9)
     *inputs, det_back = decay.map_inverse([p1, p2], [ones, ones])
-    assert np.asarray(inputs[1]) == approx(r_cos, abs=1e-8)
-    assert det * np.asarray(det_back) == approx(1.0, rel=1e-8)
+    assert_close(inputs[1], r_cos, abs=1e-8)
+    assert_close(det * det_back, 1.0, rel=1e-8)
 
 
 # --------------------------------------------------------------------------
@@ -371,12 +390,10 @@ def test_every_passing_event_is_inside_the_windows(topology, cuts):
     keep = passes(cuts, p_ext) & (det > 0)
     assert keep.sum() > 1000
     r_back, _ = windowed.map_inverse([p_ext[keep], x1[keep], x2[keep]], [])
-    r_back = np.asarray(r_back)
-    assert np.all(np.isfinite(r_back))
+    assert torch.all(torch.isfinite(r_back))
     assert r_back.min() > -1e-6 and r_back.max() < 1 + 1e-6
 
 
-
 # --------------------------------------------------------------------------
 # t channel: the first scattering, between the two beams
 # --------------------------------------------------------------------------
@@ -430,12 +447,12 @@ def test_t_channel_every_generated_point_passes():
     c = ms.Cuts(photon_cuts(AA_PIDS))
     topology = t_channel_two_body()
     _, p_free, _, _, _ = sample(ms.PhaseSpaceMapping(topology, CM_ENERGY))
-    assert np.mean(passes(c, p_free)) < 0.9
+    assert fraction(passes(c, p_free)) < 0.9
     _, p_ext, _, _, det = sample(ms.PhaseSpaceMapping(topology, CM_ENERGY, cuts=c))
-    assert np.mean(passes(c, p_ext)) > 0.999
+    assert fraction(passes(c, p_ext)) > 0.999
     photons = p_ext[:, 2:4]
-    assert np.abs(eta(photons)).max() <= ETA_PHOTON + 1e-6
-    assert np.abs(eta(photons)).max() > 0.999 * ETA_PHOTON
+    assert torch.abs(eta(photons)).max() <= ETA_PHOTON + 1e-6
+    assert torch.abs(eta(photons)).max() > 0.999 * ETA_PHOTON
 
 
 def test_t_channel_round_trip():
@@ -444,9 +461,8 @@ def test_t_channel_round_trip():
     r, p_ext, x1, x2, det = sample(mapping, n=10_000)
     keep = det > 0
     r_back, det_back = mapping.map_inverse([p_ext, x1, x2], [])
-    r_back, det_back = np.asarray(r_back), np.asarray(det_back)
-    assert r_back[keep] == approx(r[keep], abs=1e-7)
-    assert (det * det_back)[keep] == approx(1.0, rel=1e-7)
+    assert_close(r_back[keep], r[keep], abs=1e-7)
+    assert_close((det * det_back)[keep], 1.0, rel=1e-7)
 
 
 @pytest.mark.parametrize(
@@ -474,8 +490,7 @@ def test_t_channel_window_is_exact(topology, cuts):
     keep = passes(cuts, p_ext) & (det > 0)
     assert keep.sum() > 1000
     r_back, _ = windowed.map_inverse([p_ext[keep], x1[keep], x2[keep]], [])
-    r_back = np.asarray(r_back)
-    assert np.all(np.isfinite(r_back))
+    assert torch.all(torch.isfinite(r_back))
     assert r_back.min() > -1e-6 and r_back.max() < 1 + 1e-6
     assert_same_integral(topology, cuts)
 
@@ -494,9 +509,9 @@ def test_t_channel_first_step_only_bounds_what_it_peels():
     other = 2 - peeled
     cuts = ms.Cuts(photon_cuts(AAG_PIDS) + jet_cuts(AAG_PIDS, eta=None))
     _, p_ext, _, _, _ = sample(ms.PhaseSpaceMapping(topology, CM_ENERGY, cuts=cuts))
-    outside = np.abs(eta(p_ext[:, 2:5])) > ETA_PHOTON + 1e-6
-    assert outside[:, peeled].mean() < 0.03
-    assert outside[:, other].mean() > 0.2
+    outside = torch.abs(eta(p_ext[:, 2:5])) > ETA_PHOTON + 1e-6
+    assert fraction(outside[:, peeled]) < 0.03
+    assert fraction(outside[:, other]) > 0.2
 
 
 # --------------------------------------------------------------------------
@@ -515,7 +530,11 @@ def test_t_channel_first_step_only_bounds_what_it_peels():
     "1+3": [0, 2, 1, 3, 4, 5],  # double-t: {0} against {1, 2, 3}
 }
 # outgoing particles whose summed momentum has to stay within the bound
-CO_BOUNDED = {"chain": [[0], [1, 2, 3]], "2+2": [[0, 1], [2, 3]], "1+3": [[0], [1, 2, 3]]}
+CO_BOUNDED = {
+    "chain": [[0], [1, 2, 3]],
+    "2+2": [[0, 1], [2, 3]],
+    "1+3": [[0], [1, 2, 3]],
+}
 
 
 def gggg_cuts():
@@ -526,13 +545,12 @@ def gggg_cuts():
 
 
 def co_sample(mapping, n, seed):
-    rng = np.random.default_rng(seed)
-    r = rng.random((n, mapping.random_dim()))
-    d = rng.integers(0, 2, size=(n, mapping.discrete_dim())).astype(np.int32)
-    out = mapping.map_forward([r, d] if mapping.discrete_dim() else [r])
-    p_ext, x1, x2, det = (np.asarray(a) for a in out)
-    ok = np.isfinite(det) & np.all(np.isfinite(p_ext), axis=(1, 2))
-    return p_ext, x1, x2, np.where(ok, det, 0.0)
+    g = generator(seed)
+    r = rand((n, mapping.random_dim()), generator=g)
+    d = randint(0, 2, (n, mapping.discrete_dim()), generator=g, dtype=torch.int32)
+    p_ext, x1, x2, det = mapping.map_forward([r, d] if mapping.discrete_dim() else [r])
+    ok = torch.isfinite(det) & finite(p_ext)
+    return p_ext, x1, x2, torch.where(ok, det, 0.0)
 
 
 @pytest.mark.parametrize("name", list(CO_ORDERS))
@@ -546,7 +564,11 @@ def test_color_ordered_window_is_exact(name):
 
     def mapping(c):
         return ms.PhaseSpaceMapping(
-            masses, CM_ENERGY, mode=ms.PhaseSpaceMapping.color_ordered, cuts=c, color_order=order
+            masses,
+            CM_ENERGY,
+            mode=ms.PhaseSpaceMapping.color_ordered,
+            cuts=c,
+            color_order=order,
         )
 
     n = 200_000
@@ -554,13 +576,13 @@ def mapping(c):
     keep = passes(cuts, p_free) & (det_free > 0)
     assert keep.sum() > 1000
     out = mapping(cuts).map_inverse([p_free[keep], y1[keep], y2[keep]], [])
-    r_back = np.asarray(out[0])
-    assert np.all(np.isfinite(r_back))
+    r_back = out[0]
+    assert torch.all(torch.isfinite(r_back))
     assert r_back.min() > -1e-6 and r_back.max() < 1 + 1e-6
 
     def integral(p_ext, x1, x2, det):
-        w = np.where(passes(cuts, p_ext), det / (CM_ENERGY**2 * x1 * x2), 0.0)
-        return w.mean(), w.std() / math.sqrt(len(w))
+        w = torch.where(passes(cuts, p_ext), det / (CM_ENERGY**2 * x1 * x2), 0.0)
+        return mean_and_error(w)
 
     a, ea = integral(p_free, y1, y2, det_free)
     b, eb = integral(*co_sample(mapping(cuts), n, SEED + 1))
@@ -578,26 +600,27 @@ def test_color_ordered_blocks_respect_the_bounds(name):
     tau = (sqrt_s / CM_ENERGY) ** 2
     x1, x2 = math.sqrt(tau) * math.exp(y_boost), math.sqrt(tau) * math.exp(-y_boost)
     n = 100_000
-    rng = np.random.default_rng(SEED)
+    g = generator(SEED)
 
     def outside(y_max):
         mapping = ms.ColorOrderedMapping(
             order, 0.8, 0.8, [PT_JET] * 4, [], [], True, y_max
         )
-        r = rng.random((n, mapping.random_dim()))
-        d = rng.integers(0, 2, size=(n, mapping.discrete_dim())).astype(np.int32)
-        inputs = [r[:, i].copy() for i in range(r.shape[1])]
-        inputs += [d[:, j].copy() for j in range(d.shape[1])]
-        conditions = [np.full(n, sqrt_s)] + [np.zeros(n)] * 4
+        r = rand((n, mapping.random_dim()), generator=g)
+        d = randint(0, 2, (n, mapping.discrete_dim()), generator=g, dtype=torch.int32)
+        inputs = [r[:, i].clone() for i in range(r.shape[1])]
+        inputs += [d[:, j].clone() for j in range(d.shape[1])]
+        conditions = [torch.full((n,), sqrt_s)] + [torch.zeros(n)] * 4
         if y_max:
-            conditions += [np.full(n, x1), np.full(n, x2)]
+            conditions += [torch.full((n,), x1), torch.full((n,), x2)]
         out = mapping.map_forward(inputs, conditions)
-        p = np.stack([np.asarray(q) for q in out[:-1]], axis=1)
-        finite = np.all(np.isfinite(p), axis=(1, 2))
+        p = torch.stack(list(out[:-1]), dim=1)
         fractions = []
         for group in CO_BOUNDED[name]:
-            system = p[finite][:, [2 + i for i in group]].sum(axis=1)
-            fractions.append(np.mean(np.abs(rapidity(system) + y_boost) > ETA_JET + 1e-9))
+            system = p[finite(p)][:, [2 + i for i in group]].sum(dim=1)
+            fractions.append(
+                fraction(torch.abs(rapidity(system) + y_boost) > ETA_JET + 1e-9)
+            )
         return fractions
 
     for with_bound, without in zip(outside([ETA_JET] * 4), outside([])):
@@ -615,29 +638,44 @@ def test_double_t_window_alone():
     sqrt_s, y_boost = 500.0, 0.8
     tau = (sqrt_s / CM_ENERGY) ** 2
     x1, x2 = math.sqrt(tau) * math.exp(y_boost), math.sqrt(tau) * math.exp(-y_boost)
-    pa = np.tile([sqrt_s / 2, 0.0, 0.0, sqrt_s / 2], (n, 1))
-    pb = np.tile([sqrt_s / 2, 0.0, 0.0, -sqrt_s / 2], (n, 1))
-    rng = np.random.default_rng(SEED)
-    r = [rng.random(n) for _ in range(3)]
-    base = [pa, pb, np.zeros(n), np.zeros(n), np.full(n, 20.0), np.full(n, 40.0)]
+    pa = torch.tensor([sqrt_s / 2, 0.0, 0.0, sqrt_s / 2]).repeat(n, 1)
+    pb = torch.tensor([sqrt_s / 2, 0.0, 0.0, -sqrt_s / 2]).repeat(n, 1)
+    g = generator(SEED)
+    r = [rand(n, generator=g) for _ in range(3)]
+    base = [
+        pa,
+        pb,
+        torch.zeros(n),
+        torch.zeros(n),
+        torch.full((n,), 20.0),
+        torch.full((n,), 40.0),
+    ]
 
     def outside(mapping, conditions):
-        p1, p2, det = (np.asarray(a) for a in mapping.map_forward(r, conditions))
-        ok = np.isfinite(det) & (det > 0)
+        p1, p2, det = mapping.map_forward(r, conditions)
+        ok = torch.isfinite(det) & (det > 0)
         return [
-            np.mean(np.abs(rapidity(p[ok]) + y_boost) > ETA_JET + 1e-9) for p in (p1, p2)
+            fraction(torch.abs(rapidity(p[ok]) + y_boost) > ETA_JET + 1e-9)
+            for p in (p1, p2)
         ]
 
     bounded = ms.DoubleT(0.8, 0.0, 0.0, 0.8, 0.0, 0.0, True, ETA_JET, ETA_JET, 1.0)
     free = ms.DoubleT(0.8, 0.0, 0.0, 0.8, 0.0, 0.0, True)
-    assert outside(bounded, base + [np.full(n, x1), np.full(n, x2)]) == [0.0, 0.0]
+    assert outside(bounded, base + [torch.full((n,), x1), torch.full((n,), x2)]) == [
+        0.0,
+        0.0,
+    ]
     assert outside(free, base)[0] > 0.1
     *inputs, det_back = bounded.map_inverse(
-        list(bounded.map_forward(r, base + [np.full(n, x1), np.full(n, x2)])[:2]),
-        base + [np.full(n, x1), np.full(n, x2)],
+        list(
+            bounded.map_forward(r, base + [torch.full((n,), x1), torch.full((n,), x2)])[
+                :2
+            ]
+        ),
+        base + [torch.full((n,), x1), torch.full((n,), x2)],
     )
     for a, b in zip(inputs, r):
-        assert np.asarray(a) == approx(b, abs=1e-7)
+        assert_close(a, b, abs=1e-7)
 
 
 # --------------------------------------------------------------------------
@@ -665,12 +703,11 @@ def lepton_photon_cuts(pids=LA_PIDS):
 
 
 def permuted_sample(mapping, index, n, seed):
-    rng = np.random.default_rng(seed)
-    r = rng.random((n, mapping.random_dim()))
-    perm = np.full(n, index, dtype=np.int32)
-    p_ext, x1, x2, det = (np.asarray(a) for a in mapping.map_forward([r], [perm]))
-    ok = np.isfinite(det) & np.all(np.isfinite(p_ext), axis=(1, 2))
-    return r, perm, p_ext, x1, x2, np.where(ok, det, 0.0)
+    r = rand((n, mapping.random_dim()), generator=generator(seed))
+    perm = torch.full((n,), index, dtype=torch.int32)
+    p_ext, x1, x2, det = mapping.map_forward([r], [perm])
+    ok = torch.isfinite(det) & finite(p_ext)
+    return r, perm, p_ext, x1, x2, torch.where(ok, det, 0.0)
 
 
 @pytest.mark.parametrize("index", [0, 1])
@@ -688,20 +725,17 @@ def test_windows_hold_for_every_permutation(index):
     _, perm, p_free, y1, y2, det_free = permuted_sample(free, index, n, SEED)
     _, _, p_win, z1, z2, det_win = permuted_sample(windowed, index, n, SEED + 1)
 
-    assert np.abs(beam_rapidity(z1, z2)).max() <= ETA_LEPTON + 1e-9
-    assert np.mean(passes(cuts, p_win)) > np.mean(passes(cuts, p_free)) + 0.1
+    assert torch.abs(beam_rapidity(z1, z2)).max() <= ETA_LEPTON + 1e-9
+    assert fraction(passes(cuts, p_win)) > fraction(passes(cuts, p_free)) + 0.1
 
     keep = passes(cuts, p_free) & (det_free > 0)
     assert keep.sum() > 1000
-    r_back, _ = windowed.map_inverse(
-        [p_free[keep], y1[keep], y2[keep]], [perm[keep]]
-    )
-    r_back = np.asarray(r_back)
+    r_back, _ = windowed.map_inverse([p_free[keep], y1[keep], y2[keep]], [perm[keep]])
     assert r_back.min() > -1e-6 and r_back.max() < 1 + 1e-6
 
     def integral(p_ext, x1, x2, det):
-        w = np.where(passes(cuts, p_ext), det / (CM_ENERGY**2 * x1 * x2), 0.0)
-        return w.mean(), w.std() / math.sqrt(len(w))
+        w = torch.where(passes(cuts, p_ext), det / (CM_ENERGY**2 * x1 * x2), 0.0)
+        return mean_and_error(w)
 
     a, ea = integral(p_free, y1, y2, det_free)
     b, eb = integral(p_win, z1, z2, det_win)
@@ -728,9 +762,8 @@ def test_strongest_bounds_would_lose_permuted_events():
         ),
     )
     r_back, _ = too_tight.map_inverse([p_free[keep], y1[keep], y2[keep]], [])
-    r_back = np.asarray(r_back)
-    outside = ~np.all((r_back > -1e-6) & (r_back < 1 + 1e-6), axis=1)
-    assert outside.mean() > 0.01
+    outside = ~torch.all((r_back > -1e-6) & (r_back < 1 + 1e-6), dim=1)
+    assert fraction(outside) > 0.01
 
 
 # --------------------------------------------------------------------------
@@ -753,21 +786,19 @@ def test_leptonic_windows(topology):
     n = 50_000
     windowed = ms.PhaseSpaceMapping(topology, LEPTONIC_ENERGY, leptonic=True, cuts=cuts)
     free = ms.PhaseSpaceMapping(topology, LEPTONIC_ENERGY, leptonic=True)
-    rng = np.random.default_rng(SEED)
-    r = rng.random((n, windowed.random_dim()))
-    p_ext, _, _, det = (np.asarray(a) for a in windowed.map_forward([r]))
-    assert np.mean(passes(cuts, p_ext)) > 0.999
-    ones = np.ones(n)
+    r = rand((n, windowed.random_dim()), generator=generator(SEED))
+    p_ext, _, _, det = windowed.map_forward([r])
+    assert fraction(passes(cuts, p_ext)) > 0.999
+    ones = torch.ones(n)
     r_back, det_back = windowed.map_inverse([p_ext, ones, ones], [])
-    assert np.asarray(r_back) == approx(r, abs=1e-9)
-    assert det * np.asarray(det_back) == approx(1.0, rel=1e-9)
+    assert_close(r_back, r, abs=1e-9)
+    assert_close(det * det_back, 1.0, rel=1e-9)
 
-    p_free, _, _, det_free = (np.asarray(a) for a in free.map_forward([r]))
+    p_free, _, _, det_free = free.map_forward([r])
     keep = passes(cuts, p_free) & (det_free > 0)
     assert 1000 < keep.sum() < 0.9 * n
     m = int(keep.sum())
-    r_back, _ = windowed.map_inverse([p_free[keep], np.ones(m), np.ones(m)], [])
-    r_back = np.asarray(r_back)
+    r_back, _ = windowed.map_inverse([p_free[keep], torch.ones(m), torch.ones(m)], [])
     assert r_back.min() > -1e-6 and r_back.max() < 1 + 1e-6
 
 
@@ -778,10 +809,10 @@ def test_leptonic_windows(topology):
 
 def beam_conditions(n, sqrt_s=500.0, y_boost=0.8):
     tau = (sqrt_s / CM_ENERGY) ** 2
-    x1 = np.full(n, math.sqrt(tau) * math.exp(y_boost))
-    x2 = np.full(n, math.sqrt(tau) * math.exp(-y_boost))
-    pa = np.tile([sqrt_s / 2, 0.0, 0.0, sqrt_s / 2], (n, 1))
-    pb = np.tile([sqrt_s / 2, 0.0, 0.0, -sqrt_s / 2], (n, 1))
+    x1 = torch.full((n,), math.sqrt(tau) * math.exp(y_boost))
+    x2 = torch.full((n,), math.sqrt(tau) * math.exp(-y_boost))
+    pa = torch.tensor([sqrt_s / 2, 0.0, 0.0, sqrt_s / 2]).repeat(n, 1)
+    pb = torch.tensor([sqrt_s / 2, 0.0, 0.0, -sqrt_s / 2]).repeat(n, 1)
     return pa, pb, x1, x2
 
 
@@ -791,16 +822,23 @@ def test_double_t_huge_bound_is_no_bound(y_max):
     0 * inf = NaN, and the momenta are exactly those without a bound."""
     n = 20_000
     pa, pb, x1, x2 = beam_conditions(n)
-    rng = np.random.default_rng(SEED)
-    r = [rng.random(n) for _ in range(3)]
-    base = [pa, pb, np.zeros(n), np.zeros(n), np.full(n, 20.0), np.full(n, 40.0)]
+    g = generator(SEED)
+    r = [rand(n, generator=g) for _ in range(3)]
+    base = [
+        pa,
+        pb,
+        torch.zeros(n),
+        torch.zeros(n),
+        torch.full((n,), 20.0),
+        torch.full((n,), 40.0),
+    ]
     bounded = ms.DoubleT(0.8, 0.0, 0.0, 0.8, 0.0, 0.0, True, y_max, y_max, 1.0)
     free = ms.DoubleT(0.8, 0.0, 0.0, 0.8, 0.0, 0.0, True)
-    out_bounded = [np.asarray(a) for a in bounded.map_forward(r, base + [x1, x2])]
-    out_free = [np.asarray(a) for a in free.map_forward(r, base)]
+    out_bounded = bounded.map_forward(r, base + [x1, x2])
+    out_free = free.map_forward(r, base)
     for a, b in zip(out_bounded, out_free):
-        assert np.all(np.isfinite(a))
-        assert np.array_equal(a, b)
+        assert torch.all(torch.isfinite(a))
+        assert torch.equal(a, b)
 
 
 def test_loose_t_bound_leaves_the_range_untouched():
@@ -809,12 +847,12 @@ def test_loose_t_bound_leaves_the_range_untouched():
     rebuilding it through pb.k and moving its ends by rounding."""
     n = 20_000
     pa, pb, x1, x2 = beam_conditions(n)
-    rng = np.random.default_rng(SEED)
-    r = [rng.random(n) for _ in range(2)]
-    masses = [np.full(n, M_Z), np.zeros(n)]
+    g = generator(SEED)
+    r = [rand(n, generator=g) for _ in range(2)]
+    masses = [torch.full((n,), M_Z), torch.zeros(n)]
     loose = ms.TwoToTwoParticleScattering(True, 0.8, 0.0, 0.0, False, 50.0, 50.0, 1.0)
     free = ms.TwoToTwoParticleScattering(True, 0.8, 0.0, 0.0, False)
     out_loose = loose.map_forward(r + masses, [pa, pb, x1, x2])
     out_free = free.map_forward(r + masses, [pa, pb])
     for a, b in zip(out_loose, out_free):
-        assert np.array_equal(np.asarray(a), np.asarray(b))
+        assert torch.equal(a, b)
diff --git a/madspace/tests/test_cuts.py b/madspace/tests/test_cuts.py
index 952881057..2ea20ce91 100644
--- a/madspace/tests/test_cuts.py
+++ b/madspace/tests/test_cuts.py
@@ -1,8 +1,9 @@
 import math
 from dataclasses import dataclass
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand, randint
 
 import madspace as ms
 
@@ -67,11 +68,11 @@ def _color_orders(n_out):
     return orders
 
 
-def _inputs(mapping, rng, n):
-    inputs = [rng.random((n, mapping.random_dim()))]
+def _inputs(mapping, g, n):
+    inputs = [rand((n, mapping.random_dim()), generator=g)]
     dd = mapping.discrete_dim()
     if dd:
-        inputs.append(rng.integers(0, 2, size=(n, dd)).astype(np.int32))
+        inputs.append(randint(0, 2, (n, dd), generator=g, dtype=torch.int32))
     return inputs
 
 
@@ -86,26 +87,27 @@ def _cuts(pids, pt_min):
 
 
 def _feats(p):
-    pt = np.sqrt(p[..., 1] ** 2 + p[..., 2] ** 2)
-    eta = np.arcsinh(np.divide(p[..., 3], pt, out=np.zeros_like(pt), where=pt > 0))
-    phi = np.arctan2(p[..., 2], p[..., 1])
+    pt = torch.sqrt(p[..., 1] ** 2 + p[..., 2] ** 2)
+    safe_pt = torch.where(pt > 0, pt, 1.0)
+    eta = torch.arcsinh(torch.where(pt > 0, p[..., 3] / safe_pt, 0.0))
+    phi = torch.atan2(p[..., 2], p[..., 1])
     return pt, eta, phi
 
 
 def _passes(p_out, jet_idx, pt_min):
     if not jet_idx:
-        return np.ones(p_out.shape[0], dtype=bool)
+        return torch.ones(p_out.shape[0], dtype=torch.bool)
 
     pt, eta, phi = _feats(p_out)
 
-    ok = np.all(pt[:, jet_idx] >= pt_min, axis=1)
+    ok = torch.all(pt[:, jet_idx] >= pt_min, dim=1)
 
     for a in range(len(jet_idx)):
         for b in range(a + 1, len(jet_idx)):
             i, j = jet_idx[a], jet_idx[b]
-            dphi = np.abs(phi[:, i] - phi[:, j])
-            dphi = np.minimum(dphi, 2 * np.pi - dphi)
-            dR = np.sqrt((eta[:, i] - eta[:, j]) ** 2 + dphi**2)
+            dphi = torch.abs(phi[:, i] - phi[:, j])
+            dphi = torch.minimum(dphi, 2 * math.pi - dphi)
+            dR = torch.sqrt((eta[:, i] - eta[:, j]) ** 2 + dphi**2)
             ok = ok & (dR >= CUT_DELTAR_MIN)
 
     return ok
@@ -127,18 +129,18 @@ def _cut_volume(
         cuts=_cuts(pids, pt_min) if piped else None,
     )
 
-    rng = np.random.default_rng(seed)
-    p_ext, _x1, _x2, det = mapping.map_forward(_inputs(mapping, rng, n))
+    g = generator(seed)
+    p_ext, _x1, _x2, det = mapping.map_forward(_inputs(mapping, g, n))
 
-    det = np.asarray(det) * 2.0 ** mapping.discrete_dim()
-    p_out = np.asarray(p_ext)[:, 2:, :]
+    det = det * 2.0 ** mapping.discrete_dim()
+    p_out = p_ext[:, 2:, :]
 
     passes = _passes(p_out, jet_idx, pt_min)
-    finite = np.isfinite(det) & np.all(np.isfinite(p_out), axis=(1, 2))
+    finite = torch.isfinite(det) & torch.isfinite(p_out).all(dim=2).all(dim=1)
 
-    weights = np.nan_to_num(np.where(passes & finite, det, 0.0))
+    weights = torch.nan_to_num(torch.where(passes & finite, det, 0.0))
 
-    return float(weights.mean()), float(weights.std() / math.sqrt(n))
+    return weights.mean().item(), weights.std(correction=0).item() / math.sqrt(n)
 
 
 def _configs(n_out):
diff --git a/madspace/tests/test_decay_topology.py b/madspace/tests/test_decay_topology.py
index a58fc9c9b..6ed5ec0c8 100644
--- a/madspace/tests/test_decay_topology.py
+++ b/madspace/tests/test_decay_topology.py
@@ -13,8 +13,9 @@
 
 import math
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand
 from pytest import approx
 
 import madspace as ms
@@ -25,7 +26,13 @@
 M_W = 80.379
 W_W = 2.085
 
-rng = np.random.default_rng(20260811)
+rng = generator(20260812)
+
+
+def assert_close(actual, expected, rtol=1e-6, atol=1e-12):
+    """pytest.approx for tensors; expected may be a number."""
+    expected = torch.as_tensor(expected, dtype=actual.dtype).expand_as(actual)
+    torch.testing.assert_close(actual, expected, rtol=rtol, atol=atol)
 
 
 def massless_volume(n, mass):
@@ -54,20 +61,21 @@ def decay_mapping(outgoing_masses, propagators, vertices, mass=M_TOP):
 
 
 def sample(mapping):
-    r = rng.random((BATCH_SIZE, mapping.random_dim()))
+    r = rand((BATCH_SIZE, mapping.random_dim()), generator=rng)
     p_ext, x1, x2, det = mapping.map_forward([r], [])
     return r, p_ext, x1, x2, det
 
 
 def volume(det):
     """Monte-Carlo estimate of the volume and its standard error."""
-    return np.mean(det), np.std(det) / math.sqrt(len(det))
+    return det.mean().item(), det.std(correction=0).item() / math.sqrt(len(det))
 
 
 # --------------------------------------------------------------------------
 # Topology
 # --------------------------------------------------------------------------
 
+
 def test_single_incoming_has_no_t_channel():
     diagram = ms.Diagram([M_TOP], [0.0, 0.0], [], [["i0", "o0", "o1"]])
     topology = ms.Topology(diagram)
@@ -85,13 +93,14 @@ def test_three_or_more_incoming_rejected():
 # Kinematics
 # --------------------------------------------------------------------------
 
+
 def test_decaying_particle_is_at_rest():
     mapping = decay_mapping([0.0, 0.0], [], [["i0", "o0", "o1"]])
     _, p_ext, _, _, _ = sample(mapping)
     # One incoming + two outgoing, not the 2 + n of a collision.
     assert p_ext.shape == (BATCH_SIZE, 3, 4)
-    assert p_ext[:, 0, 0] == approx(M_TOP)
-    assert p_ext[:, 0, 1:] == approx(0.0, abs=1e-10)
+    assert_close(p_ext[:, 0, 0], M_TOP)
+    assert_close(p_ext[:, 0, 1:], 0.0, rtol=0.0, atol=1e-10)
 
 
 def test_momentum_conservation_and_masses():
@@ -102,23 +111,23 @@ def test_momentum_conservation_and_masses():
     )
     _, p_ext, _, _, _ = sample(mapping)
     assert p_ext.shape == (BATCH_SIZE, 4, 4)
-    assert np.sum(p_ext[:, 1:], axis=1) == approx(p_ext[:, 0], abs=1e-8)
+    assert_close(torch.sum(p_ext[:, 1:], dim=1), p_ext[:, 0], rtol=0.0, atol=1e-8)
 
     def mass(p):
-        m2 = p[:, 0] ** 2 - np.sum(p[:, 1:] ** 2, axis=1)
-        return np.sqrt(np.abs(m2))
+        m2 = p[:, 0] ** 2 - torch.sum(p[:, 1:] ** 2, dim=1)
+        return torch.sqrt(torch.abs(m2))
 
-    assert mass(p_ext[:, 0]) == approx(M_TOP, rel=1e-10)
-    assert mass(p_ext[:, 1]) == approx(M_B, rel=1e-6)
-    assert mass(p_ext[:, 2]) == approx(0.0, abs=1e-5)
-    assert mass(p_ext[:, 3]) == approx(0.0, abs=1e-5)
+    assert_close(mass(p_ext[:, 0]), M_TOP, rtol=1e-10)
+    assert_close(mass(p_ext[:, 1]), M_B, rtol=1e-6)
+    assert_close(mass(p_ext[:, 2]), 0.0, rtol=0.0, atol=1e-5)
+    assert_close(mass(p_ext[:, 3]), 0.0, rtol=0.0, atol=1e-5)
 
 
 def test_no_beam_momentum_fractions():
     mapping = decay_mapping([0.0, 0.0], [], [["i0", "o0", "o1"]])
     _, _, x1, x2, _ = sample(mapping)
-    assert x1 == approx(1.0)
-    assert x2 == approx(1.0)
+    assert_close(x1, 1.0)
+    assert_close(x2, 1.0)
 
 
 def test_random_dim_is_3n_minus_4():
@@ -126,9 +135,11 @@ def test_random_dim_is_3n_minus_4():
         mapping = decay_mapping(
             [0.0] * n,
             [] if n == 2 else [ms.Propagator(mass=M_W, width=W_W, pdg_id=24)],
-            [["i0", "o0", "o1"]]
-            if n == 2
-            else [["i0", "o0", "p0"], ["p0", "o1", "o2"]],
+            (
+                [["i0", "o0", "o1"]]
+                if n == 2
+                else [["i0", "o0", "p0"], ["p0", "o1", "o2"]]
+            ),
         )
         assert mapping.random_dim() == 3 * n - 4
         assert mapping.particle_count() == n + 1
@@ -138,6 +149,7 @@ def test_random_dim_is_3n_minus_4():
 # Phase-space volume
 # --------------------------------------------------------------------------
 
+
 def test_two_body_massless_volume():
     mapping = decay_mapping([0.0, 0.0], [], [["i0", "o0", "o1"]])
     _, _, _, _, det = sample(mapping)
@@ -185,6 +197,7 @@ def test_three_body_massless_volume(prop_mass, prop_width, rel_tol):
 # Invertibility
 # --------------------------------------------------------------------------
 
+
 @pytest.mark.parametrize(
     "outgoing,propagators,vertices",
     [
@@ -202,5 +215,5 @@ def test_forward_inverse_round_trip(outgoing, propagators, vertices):
     mapping = decay_mapping(outgoing, propagators, vertices)
     r, p_ext, x1, x2, det = sample(mapping)
     r_back, det_back = mapping.map_inverse([p_ext, x1, x2], [])
-    assert r_back == approx(r, abs=1e-8)
-    assert det * det_back == approx(1.0, rel=1e-8)
+    assert_close(r_back, r, rtol=0.0, atol=1e-8)
+    assert_close(det * det_back, 1.0, rtol=1e-8)
diff --git a/madspace/tests/test_discrete_sampler.py b/madspace/tests/test_discrete_sampler.py
index 4f5b8c73e..d8b29508a 100644
--- a/madspace/tests/test_discrete_sampler.py
+++ b/madspace/tests/test_discrete_sampler.py
@@ -11,7 +11,8 @@
 need to special-case a "degenerate" result.
 """
 
-import numpy as np
+import torch
+from madspace_testing import default_context
 
 import madspace as ms
 
@@ -20,15 +21,15 @@
 # row 0: a healthy prior
 # row 1: every option's density exactly zero
 # row 2: a single non-zero option
-PRIOR = np.array(
+PRIOR = torch.tensor(
     [
         [0.3, 0.2, 0.4, 0.1],
         [0.0, 0.0, 0.0, 0.0],
         [0.0, 0.0, 0.0, 0.7],
     ]
 )
-RANDOM = np.array([0.42, 0.42, 0.42])
-UNUSED_INDEX = np.zeros(len(RANDOM), dtype=np.int32)
+RANDOM = torch.tensor([0.42, 0.42, 0.42])
+UNUSED_INDEX = torch.zeros(len(RANDOM), dtype=torch.int32)
 ALL_ZERO_ROW = 1
 
 
@@ -40,7 +41,7 @@ def sampler_and_context():
     # own name for the probability global
     global _prefix_count
     _prefix_count += 1
-    ctx = ms.default_context()
+    ctx = default_context()
     sampler = ms.DiscreteSampler([N_OPTIONS], f"test_zero_prob{_prefix_count}_", [0])
     sampler.initialize_globals(ctx)
     return sampler, ctx
@@ -63,7 +64,7 @@ def run(builder_fn, output_names, inputs, output_types=None):
     )
     builder_fn(fb, sampler)
     runtime = ms.FunctionRuntime(fb.function(), ctx)
-    return dict(zip(output_names, (ms.Tensor.numpy(t) for t in runtime.call(inputs))))
+    return dict(zip(output_names, (torch.from_dlpack(t) for t in runtime.call(inputs))))
 
 
 def test_forward_det_is_always_at_least_one():
@@ -85,7 +86,7 @@ def build(fb, sampler):
         [RANDOM, PRIOR, UNUSED_INDEX],
         output_types={"index": ms.batch_int},
     )
-    assert np.all(out["det"] >= 1.0)
+    assert torch.all(out["det"] >= 1.0)
     assert out["det"][ALL_ZERO_ROW] == 1.0
     assert out["index"][ALL_ZERO_ROW] == 0
 
@@ -95,9 +96,9 @@ def test_forward_selects_the_expected_bucket():
     cumulative-probability bucket contains it. For [0.3, 0.2, 0.4, 0.1] the
     bucket edges are [0.3, 0.5, 0.9, 1.0].
     """
-    prior = np.tile([0.3, 0.2, 0.4, 0.1], (5, 1))
-    random = np.array([0.0, 0.29, 0.31, 0.89, 0.99])
-    expected = np.array([0, 0, 1, 2, 3])
+    prior = torch.tensor([0.3, 0.2, 0.4, 0.1]).repeat(5, 1)
+    random = torch.tensor([0.0, 0.29, 0.31, 0.89, 0.99])
+    expected = torch.tensor([0, 0, 1, 2, 3], dtype=torch.int32)
 
     def build(fb, sampler):
         fb.output(0, sampler.build_forward(fb, [fb.input(0)], [fb.input(1)])[0])
@@ -108,7 +109,7 @@ def build(fb, sampler):
         [random, prior, UNUSED_INDEX[: len(random)]],
         output_types={"index": ms.batch_int},
     )
-    assert np.array_equal(out["index"], expected)
+    assert torch.equal(out["index"], expected)
 
 
 def test_forward_never_selects_a_zero_probability_option():
@@ -117,8 +118,8 @@ def test_forward_never_selects_a_zero_probability_option():
     directly, rather than relying on r never landing exactly on their
     zero-width buckets.
     """
-    prior = np.tile([0.0, 0.0, 0.0, 0.7], (5, 1))
-    random = np.array([0.0, 0.25, 0.5, 0.75, 0.999999])
+    prior = torch.tensor([0.0, 0.0, 0.0, 0.7]).repeat(5, 1)
+    random = torch.tensor([0.0, 0.25, 0.5, 0.75, 0.999999])
 
     def build(fb, sampler):
         fb.output(0, sampler.build_forward(fb, [fb.input(0)], [fb.input(1)])[0])
@@ -129,7 +130,7 @@ def build(fb, sampler):
         [random, prior, UNUSED_INDEX[: len(random)]],
         output_types={"index": ms.batch_int},
     )
-    assert np.all(out["index"] == 3)
+    assert torch.all(out["index"] == 3)
 
 
 def test_forward_and_inverse_dets_are_reciprocal():
@@ -147,7 +148,8 @@ def build(fb, sampler):
         fb.output(1, inverse["det"])
 
     out = run(build, ["det_fwd", "det_inv"], [RANDOM, PRIOR, UNUSED_INDEX])
-    assert np.allclose(out["det_fwd"] * out["det_inv"], 1.0)
+    product = out["det_fwd"] * out["det_inv"]
+    assert torch.allclose(product, torch.ones_like(product))
 
 
 def test_inverse_all_zero_prior_returns_the_canonical_draw():
@@ -156,14 +158,14 @@ def test_inverse_all_zero_prior_returns_the_canonical_draw():
     det = 1 no matter which index is asked for, matching the forward
     sampler's index-0/probability-1 convention for the same row.
     """
-    prior = np.zeros((3, N_OPTIONS))
-    index_in = np.array([0, 1, 3], dtype=np.int32)
+    prior = torch.zeros((3, N_OPTIONS))
+    index_in = torch.tensor([0, 1, 3], dtype=torch.int32)
 
     def build(fb, sampler):
         inverse = sampler.build_inverse(fb, [fb.input(2)], [fb.input(1)])
         fb.output(0, inverse[0])
         fb.output(1, inverse["det"])
 
-    out = run(build, ["r", "det"], [np.zeros(3), prior, index_in])
-    assert np.all(out["r"] == 0.5)
-    assert np.all(out["det"] == 1.0)
+    out = run(build, ["r", "det"], [torch.zeros(3), prior, index_in])
+    assert torch.all(out["r"] == 0.5)
+    assert torch.all(out["det"] == 1.0)
diff --git a/madspace/tests/test_double_t.py b/madspace/tests/test_double_t.py
index 4abbd26a5..35de0fbbe 100644
--- a/madspace/tests/test_double_t.py
+++ b/madspace/tests/test_double_t.py
@@ -1,8 +1,9 @@
+import math
 from collections import namedtuple
 
-import numpy as np
 import pytest
-from pytest import approx
+import torch
+from madspace_testing import generator, rand, randn
 
 import madspace as ms
 
@@ -12,11 +13,38 @@
 
 
 def inv_mass_sq(p):
-    return p[..., 0] ** 2 - np.sum(p[..., 1:] ** 2, axis=-1)
+    return p[..., 0] ** 2 - torch.sum(p[..., 1:] ** 2, dim=-1)
 
 
 def mass(p):
-    return np.sqrt(np.maximum(inv_mass_sq(p), 0.0))
+    return torch.sqrt(torch.clamp(inv_mass_sq(p), min=0.0))
+
+
+def uniform(g, a, b, n):
+    return a + (b - a) * rand(n, generator=g)
+
+
+def beams(pz):
+    """Massless beams (pz, 0, 0, +pz) and (pz, 0, 0, -pz)."""
+    zeros = torch.zeros_like(pz)
+    return (
+        torch.stack([pz, zeros, zeros, +pz], dim=1),
+        torch.stack([pz, zeros, zeros, -pz], dim=1),
+    )
+
+
+def lab_beams(g):
+    """Massless beams with momenta smeared around (0, 0, +-4000)."""
+    pa3 = torch.tensor([[0, 0, 4000.0]]) + 500.0 * randn((N, 3), generator=g)
+    ea = torch.sqrt(torch.sum(pa3**2, dim=1))  # massless
+    pb3 = torch.tensor([[0, 0, -4000.0]]) + 500.0 * randn((N, 3), generator=g)
+    eb = torch.sqrt(torch.sum(pb3**2, dim=1))
+    return torch.cat([ea[:, None], pa3], dim=1), torch.cat([eb[:, None], pb3], dim=1)
+
+
+def assert_approx(actual, expected, rel=1e-6, abs=1e-12, msg=None):
+    """pytest.approx-like elementwise comparison of tensors."""
+    torch.testing.assert_close(actual, expected, rtol=rel, atol=abs, msg=msg)
 
 
 # ----------------------------
@@ -41,7 +69,7 @@ def mass(p):
 
 @pytest.fixture
 def rng():
-    return np.random.default_rng(1234)
+    return generator(1234)
 
 
 @pytest.fixture(
@@ -55,39 +83,31 @@ def input_points(rng, request):
     """
     com = request.param
 
-    pz = rng.uniform(1000.0, 4000.0, N)
-    e1 = pz  # massless beams in COM
-    e2 = pz
-    pa = np.stack([e1, np.zeros(N), np.zeros(N), +pz], axis=1)
-    pb = np.stack([e2, np.zeros(N), np.zeros(N), -pz], axis=1)
+    pz = uniform(rng, 1000.0, 4000.0, N)
+    pa, pb = beams(pz)  # massless beams in COM
 
     if not com:
         # boost to a random lab frame (the same recipe used in the 2->3 test)
-        pa3 = np.array([[0, 0, 4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-        ea = np.sqrt(np.sum(pa3**2, axis=1))  # massless
-        pa = np.concatenate([ea[:, None], pa3], axis=1)
-        pb3 = np.array([[0, 0, -4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-        eb = np.sqrt(np.sum(pb3**2, axis=1))
-        pb = np.concatenate([eb[:, None], pb3], axis=1)
+        pa, pb = lab_beams(rng)
 
     # m1: the fixed single-particle mass
     # mir_min: lower bound on recoil mass (e.g. sum of constituent masses).
     # Pick mir_min well below the sqrt(s) so there\'s room for the recoil.
-    m1 = rng.uniform(1.0, 50.0, N)
-    mir_min = rng.uniform(50.0, 200.0, N)
+    m1 = uniform(rng, 1.0, 50.0, N)
+    mir_min = uniform(rng, 50.0, 200.0, N)
 
-    r_phi = rng.random(N)
-    r_t1 = rng.random(N)
-    r_t2 = rng.random(N)
+    r_phi = rand(N, generator=rng)
+    r_t1 = rand(N, generator=rng)
+    r_t2 = rand(N, generator=rng)
 
     return InputPoint(r_phi, r_t1, r_t2, m1, mir_min, pa, pb)
 
 
 @pytest.fixture(
     params=[
-        (np.full(N, 5.0), np.full(N, 80.0)),
-        (np.full(N, 50.0), np.full(N, 150.0)),
-        (np.full(N, 100.0), np.full(N, 250.0)),
+        (5.0, 80.0),
+        (50.0, 150.0),
+        (100.0, 250.0),
     ],
     ids=[
         "light single, light recoil min",
@@ -97,17 +117,13 @@ def input_points(rng, request):
 )
 def fixed_input_points(rng, request):
     """Fixed (m1, mir_min) used by the phase-space-volume test."""
-    m1, mir_min = request.param
+    m1, mir_min = (torch.full((N,), m) for m in request.param)
 
-    pz = np.full(N, 3000.0)
-    e1 = pz
-    e2 = pz
-    pa = np.stack([e1, np.zeros(N), np.zeros(N), +pz], axis=1)
-    pb = np.stack([e2, np.zeros(N), np.zeros(N), -pz], axis=1)
+    pa, pb = beams(torch.full((N,), 3000.0))
 
-    r_phi = rng.random(N)
-    r_t1 = rng.random(N)
-    r_t2 = rng.random(N)
+    r_phi = rand(N, generator=rng)
+    r_t1 = rand(N, generator=rng)
+    r_t2 = rand(N, generator=rng)
 
     return InputPoint(r_phi, r_t1, r_t2, m1, mir_min, pa, pb)
 
@@ -133,7 +149,7 @@ def test_momentum_conservation(input_points):
     ]
 
     p1, p2, det = mapping.map_forward(inputs, conditions)
-    assert p1 + p2 == approx(input_points.pa + input_points.pb)
+    assert_approx(p1 + p2, input_points.pa + input_points.pb)
 
 
 def test_on_shell_mass(input_points):
@@ -152,7 +168,7 @@ def test_on_shell_mass(input_points):
     ]
 
     p1, p2, det = mapping.map_forward(inputs, conditions)
-    assert mass(p1) == approx(input_points.m1)
+    assert_approx(mass(p1), input_points.m1)
 
 
 def test_recoil_mass_minimum(input_points):
@@ -175,7 +191,7 @@ def test_recoil_mass_minimum(input_points):
     p1, p2, det = mapping.map_forward(inputs, conditions)
     m_ir = mass(p2)
     # allow a small numerical tolerance below mir_min
-    assert np.all(m_ir >= input_points.mir_min - 1e-6 * input_points.mir_min)
+    assert torch.all(m_ir >= input_points.mir_min - 1e-6 * input_points.mir_min)
 
 
 def test_inverse(input_points):
@@ -196,10 +212,10 @@ def test_inverse(input_points):
     p1, p2, det = mapping.map_forward(inputs, conditions)
     *inv_inputs, inv_det = mapping.map_inverse([p1, p2], conditions)
 
-    assert inv_det == approx(1 / det, rel=1e-5)
+    assert_approx(inv_det, 1 / det, rel=1e-5)
 
     for i, (inp, inv_inp) in enumerate(zip(inputs, inv_inputs)):
-        assert inp == approx(inv_inp), f"mismatch in input index {i}"
+        assert_approx(inp, inv_inp, msg=f"mismatch in input index {i}")
 
 
 def test_phase_space_volume_analytic(fixed_input_points):
@@ -230,9 +246,9 @@ def test_phase_space_volume_analytic(fixed_input_points):
     # All entries of fixed_input_points have the same s, m1, mir_min by construction.
     pa = fixed_input_points.pa[0]
     pb = fixed_input_points.pb[0]
-    s = (pa[0] + pb[0]) ** 2 - np.sum((pa[1:] + pb[1:]) ** 2)
-    m1_val = fixed_input_points.m1[0]
-    mir_min_val = fixed_input_points.mir_min[0]
+    s = float((pa[0] + pb[0]) ** 2 - torch.sum((pa[1:] + pb[1:]) ** 2))
+    m1_val = float(fixed_input_points.m1[0])
+    mir_min_val = float(fixed_input_points.mir_min[0])
 
     def lam(s, m1_sq, m2_sq):
         return (s - m1_sq - m2_sq) ** 2 - 4 * m1_sq * m2_sq
@@ -241,16 +257,16 @@ def integrand(mir_sq):
         l = lam(s, m1_val**2, mir_sq)
         if l <= 0:
             return 0.0
-        return np.pi / (2 * s) * np.sqrt(l)
+        return math.pi / (2 * s) * math.sqrt(l)
 
     analytic_volume, _ = quad(
         integrand,
         mir_min_val**2,
-        (np.sqrt(s) - m1_val) ** 2,
+        (math.sqrt(s) - m1_val) ** 2,
     )
 
-    mc_volume = np.mean(det)
-    mc_err = np.std(det) / np.sqrt(N)
+    mc_volume = det.mean().item()
+    mc_err = det.std(correction=0).item() / math.sqrt(N)
 
     # 5-sigma agreement
     diff = abs(mc_volume - analytic_volume)
@@ -271,19 +287,14 @@ def test_inverse_lab_specifically(rng):
     mapping = ms.DoubleT()
 
     # boosted lab frame
-    pz = rng.uniform(1000.0, 4000.0, N)
-    pa3 = np.array([[0, 0, 4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-    ea = np.sqrt(np.sum(pa3**2, axis=1))
-    pa = np.concatenate([ea[:, None], pa3], axis=1)
-    pb3 = np.array([[0, 0, -4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-    eb = np.sqrt(np.sum(pb3**2, axis=1))
-    pb = np.concatenate([eb[:, None], pb3], axis=1)
-
-    m1 = np.full(N, 30.0)
-    mir_min = np.full(N, 100.0)
-    r_phi = rng.random(N)
-    r_t1 = rng.random(N)
-    r_t2 = rng.random(N)
+    pz = uniform(rng, 1000.0, 4000.0, N)
+    pa, pb = lab_beams(rng)
+
+    m1 = torch.full((N,), 30.0)
+    mir_min = torch.full((N,), 100.0)
+    r_phi = rand(N, generator=rng)
+    r_t1 = rand(N, generator=rng)
+    r_t2 = rand(N, generator=rng)
 
     inputs = [r_phi, r_t1, r_t2]
     conditions = [pa, pb, m1, mir_min]
@@ -293,8 +304,8 @@ def test_inverse_lab_specifically(rng):
 
     names = ["r_phi", "r_t1", "r_t2"]
     for name, inp, inv_inp in zip(names, inputs, inv_inputs):
-        assert inp == approx(inv_inp), f"LAB mismatch in {name}"
-    assert inv_det == approx(1 / det, rel=1e-5)
+        assert_approx(inp, inv_inp, msg=f"LAB mismatch in {name}")
+    assert_approx(inv_det, 1 / det, rel=1e-5)
 
 
 # ----------------------------
@@ -317,16 +328,14 @@ def massless_input_points(rng, request):
     condition (Option B fix), this should round-trip cleanly."""
     _, m1_val = request.param
 
-    pz = np.full(N, 6500.0)  # 13 TeV
-    pa = np.stack([pz, np.zeros(N), np.zeros(N), +pz], axis=1)
-    pb = np.stack([pz, np.zeros(N), np.zeros(N), -pz], axis=1)
+    pa, pb = beams(torch.full((N,), 6500.0))  # 13 TeV
 
-    m1 = np.full(N, m1_val)
-    mir_min = np.full(N, 50.0)
+    m1 = torch.full((N,), m1_val)
+    mir_min = torch.full((N,), 50.0)
 
-    r_phi = rng.random(N)
-    r_t1 = rng.random(N)
-    r_t2 = rng.random(N)
+    r_phi = rand(N, generator=rng)
+    r_t1 = rand(N, generator=rng)
+    r_t2 = rand(N, generator=rng)
 
     return InputPoint(r_phi, r_t1, r_t2, m1, mir_min, pa, pb)
 
@@ -351,14 +360,14 @@ def test_inverse_massless(massless_input_points):
     *inv_inputs, inv_det = mapping.map_inverse([p1, p2], conditions)
 
     # det round-trip
-    assert inv_det == approx(1 / det, rel=1e-5)
+    assert_approx(inv_det, 1 / det, rel=1e-5)
 
     # random round-trip; r_t2 used to be the most sensitive (its bounds
     # depend on t1_abs which depended on the lossy m1^2 recompute).
     for i, (inp, inv_inp) in enumerate(zip(inputs, inv_inputs)):
-        assert inp == approx(
-            inv_inp, rel=1e-5, abs=1e-8
-        ), f"mismatch in input index {i}"
+        assert_approx(
+            inp, inv_inp, rel=1e-5, abs=1e-8, msg=f"mismatch in input index {i}"
+        )
 
 
 def test_phase_space_volume_massless(rng):
@@ -369,15 +378,15 @@ def test_phase_space_volume_massless(rng):
     mir_min_val = 80.0
     s_val = 13000.0**2
 
-    pa = np.tile([6500.0, 0, 0, 6500.0], (N, 1))
-    pb = np.tile([6500.0, 0, 0, -6500.0], (N, 1))
+    pa = torch.tensor([6500.0, 0, 0, 6500.0]).repeat(N, 1)
+    pb = torch.tensor([6500.0, 0, 0, -6500.0]).repeat(N, 1)
 
-    m1 = np.full(N, m1_val)
-    mir_min = np.full(N, mir_min_val)
+    m1 = torch.full((N,), m1_val)
+    mir_min = torch.full((N,), mir_min_val)
 
-    r_phi = rng.random(N)
-    r_t1 = rng.random(N)
-    r_t2 = rng.random(N)
+    r_phi = rand(N, generator=rng)
+    r_t1 = rand(N, generator=rng)
+    r_t2 = rand(N, generator=rng)
 
     mapping = ms.DoubleT()
     p1, p2, det = mapping.map_forward(
@@ -392,16 +401,16 @@ def integrand(mir_sq):
         l = lam(s_val, m1_val**2, mir_sq)
         if l <= 0:
             return 0.0
-        return np.pi / (2 * s_val) * np.sqrt(l)
+        return math.pi / (2 * s_val) * math.sqrt(l)
 
     analytic_volume, _ = quad(
         integrand,
         mir_min_val**2,
-        (np.sqrt(s_val) - m1_val) ** 2,
+        (math.sqrt(s_val) - m1_val) ** 2,
     )
 
-    mc_volume = np.mean(det)
-    mc_err = np.std(det) / np.sqrt(N)
+    mc_volume = det.mean().item()
+    mc_err = det.std(correction=0).item() / math.sqrt(N)
 
     diff = abs(mc_volume - analytic_volume)
     assert diff < 5 * mc_err, (
diff --git a/madspace/tests/test_event_histograms.py b/madspace/tests/test_event_histograms.py
index adaf490ed..65873dd71 100644
--- a/madspace/tests/test_event_histograms.py
+++ b/madspace/tests/test_event_histograms.py
@@ -11,10 +11,14 @@
 import json
 
 import numpy as np
+import pytest
 from pytest import approx
 
 import madspace as ms
 
+# host-side: EventHistograms fills from python lists on an explicit CPU context
+pytestmark = pytest.mark.cpu_only
+
 
 def random_momenta(n, seed=5):
     """2 -> 2 massless momenta (E px py pz), incoming along z."""
@@ -40,6 +44,7 @@ def boost(e, pz):
         out.append(p_in + [[ea, px, py, pza], [eb, -px, -py, pzb]])
     return out
 
+
 def test_weight_histogram():
     """A from_weight histogram bins the event weight over the reference one,
     so an unweighted sample is a spike at 1 whatever the cross section."""
@@ -47,11 +52,15 @@ def test_weight_histogram():
     momenta = random_momenta(n, seed=11)
     ctx = ms.Context(device=ms.cpu_device(), thread_count=1)
     pids = [2, 2, 2, 2]
-    obs = ms.ObservableValues([
-        ms.Observable(pids, observable="sqrt_s", select_pids=[], name="sqrt_s"),
-    ])
-    specs = [ms.EventHistogramSpec("sqrt_s", 0.0, 2000.0, 8),
-             ms.EventHistogramSpec("weight", 0.0, 5.0, 10, from_weight=True)]
+    obs = ms.ObservableValues(
+        [
+            ms.Observable(pids, observable="sqrt_s", select_pids=[], name="sqrt_s"),
+        ]
+    )
+    specs = [
+        ms.EventHistogramSpec("sqrt_s", 0.0, 2000.0, 8),
+        ms.EventHistogramSpec("weight", 0.0, 5.0, 10, from_weight=True),
+    ]
     # half the events carry twice the reference weight
     weight = np.full(n, 3.0)
     weight[::2] = 6.0
@@ -84,4 +93,4 @@ def test_weight_histogram_alone():
     data = json.loads(hists.to_json())
     values = data[0]["bin_values"]
     assert sum(values) == approx(2.0)
-    assert values[3] == approx(2.0)   # the bin holding w/<w> = 1
+    assert values[3] == approx(2.0)  # the bin holding w/<w> = 1
diff --git a/madspace/tests/test_flow.py b/madspace/tests/test_flow.py
index d9b302a7a..f7ee2acbc 100644
--- a/madspace/tests/test_flow.py
+++ b/madspace/tests/test_flow.py
@@ -2,14 +2,12 @@
 import torch
 import torch.nn as nn
 import torch.nn.functional as F
+from madspace_testing import default_context, generator, rand, randn
 from pytest import approx
 
 import madspace as ms
 from madspace.torch import FunctionModule
 
-torch.set_default_dtype(torch.float64)
-torch.manual_seed(3210)
-
 
 @pytest.mark.parametrize("inverse", [False, True], ids=["forward", "inverse"])
 def test_block_gradient(inverse):
@@ -42,9 +40,10 @@ def test_block_gradient(inverse):
     fb.output(0, out)
     fb.output(1, fb.reduce_product(det))
     func = fb.function()
-    module = FunctionModule(func)
+    module = FunctionModule(func, default_context())
 
     n_points = 200
-    r_in = torch.rand((n_points, n_dims), requires_grad=True)
-    r_cond = torch.randn((n_points, n_cond), requires_grad=True)
+    g = generator(3210)
+    r_in = rand((n_points, n_dims), generator=g, requires_grad=True)
+    r_cond = randn((n_points, n_cond), generator=g, requires_grad=True)
     torch.autograd.gradcheck(module, [r_in, r_cond], raise_exception=True, rtol=1e-4)
diff --git a/madspace/tests/test_histogram_op.py b/madspace/tests/test_histogram_op.py
index 3cc026c3c..f186cd91f 100644
--- a/madspace/tests/test_histogram_op.py
+++ b/madspace/tests/test_histogram_op.py
@@ -2,45 +2,87 @@
 
 Its outputs hold n_bins + 2 values, in that order: underflow, the interior bins, overflow.
 The GPU runtime reduces the events to these n_bins + 2 keys (and used to do so into a
-scratch buffer of n_bins); .github/gpu_ci/madspace_checks.py compares it with this CPU
-result on a cluster GPU.
+scratch buffer of n_bins, summing the squared weights of the wrong events in each bin);
+the GPU CI runs these tests on a GPU via MADSPACE_TEST_DEVICE, where the large histogram
+is compared with the CPU runtime.
 """
 
-import numpy as np
+import torch
+from madspace_testing import default_context, generator, rand, randn
 
 import madspace as ms
 
 N_BINS = 5
 
 
-def run_histogram(x, w, lo, hi):
-    out_type = ms.Type(ms.DataType.float, ms.BatchSize.one, [N_BINS + 2])
+def histogram_function(n_bins):
+    out_type = ms.Type(ms.DataType.float, ms.BatchSize.one, [n_bins + 2])
     fb = ms.FunctionBuilder(
         ms.NamedTypes([(name, ms.batch_float) for name in ("x", "w", "lo", "hi")]),
         ms.NamedTypes([("values", out_type), ("square_values", out_type)]),
     )
-    values, square_values = fb.histogram(*(fb.input(i) for i in range(4)), N_BINS)
+    values, square_values = fb.histogram(*(fb.input(i) for i in range(4)), n_bins)
     fb.output(0, values)
     fb.output(1, square_values)
-    runtime = ms.FunctionRuntime(fb.function(), ms.default_context())
-    return [np.asarray(ms.Tensor.numpy(t))[0] for t in runtime.call([x, w, lo, hi])]
+    return fb.function()
+
+
+def run_histogram(x, w, lo, hi):
+    runtime = ms.FunctionRuntime(histogram_function(N_BINS), default_context())
+    return [torch.from_dlpack(t)[0] for t in runtime.call([x, w, lo, hi])]
 
 
 def test_histogram_underflow_bins_overflow():
-    x = np.array([-1.0, 0.1, 0.5, 0.95, 1.5, 2.0, 0.3])
-    w = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0])
+    x = torch.tensor([-1.0, 0.1, 0.5, 0.95, 1.5, 2.0, 0.3])
+    w = torch.tensor([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0])
     n = len(x)
-    values, square_values = run_histogram(x, w, np.zeros(n), np.ones(n))
+    values, square_values = run_histogram(x, w, torch.zeros(n), torch.ones(n))
     # underflow: -1; bins of width 0.2: 0.1, 0.3, 0.5, -, 0.95; overflow: 1.5 and 2
-    expected = np.array([1.0, 2.0, 7.0, 3.0, 0.0, 4.0, 5.0 + 6.0])
-    np.testing.assert_allclose(values, expected)
-    np.testing.assert_allclose(square_values, [1.0, 4.0, 49.0, 9.0, 0.0, 16.0, 25.0 + 36.0])
+    expected = torch.tensor([1.0, 2.0, 7.0, 3.0, 0.0, 4.0, 5.0 + 6.0])
+    torch.testing.assert_close(values, expected, rtol=1e-7, atol=0)
+    torch.testing.assert_close(
+        square_values,
+        torch.tensor([1.0, 4.0, 49.0, 9.0, 0.0, 16.0, 25.0 + 36.0]),
+        rtol=1e-7,
+        atol=0,
+    )
 
 
 def test_histogram_range_per_event():
     """min and max are per-event inputs: the same value falls in different bins."""
-    x = np.full(3, 0.5)
-    w = np.ones(3)
-    values, _ = run_histogram(x, w, np.array([0.0, 0.4, 1.0]), np.array([1.0, 1.4, 2.0]))
+    x = torch.full((3,), 0.5)
+    w = torch.ones(3)
+    values, _ = run_histogram(
+        x, w, torch.tensor([0.0, 0.4, 1.0]), torch.tensor([1.0, 1.4, 2.0])
+    )
     # 0.5 in [0, 1): bin 2 (index 3); in [0.4, 1.4): bin 0 (index 1); below [1, 2): underflow
-    np.testing.assert_allclose(values, [1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0])
+    torch.testing.assert_close(
+        values, torch.tensor([1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0]), rtol=1e-7, atol=0
+    )
+
+
+def test_large_histogram_matches_cpu_runtime(n=100000, n_bins=40):
+    """Many events in many bins on the test device, against the CPU runtime."""
+    function = histogram_function(n_bins)
+    g = generator(17)
+    inputs = [
+        -0.5 + 2.0 * rand(n, generator=g),
+        randn(n, generator=g),
+        torch.zeros(n),
+        torch.ones(n),
+    ]
+    device = [
+        torch.from_dlpack(t).cpu()
+        for t in ms.FunctionRuntime(function, default_context()).call(inputs)
+    ]
+    cpu = [
+        torch.from_dlpack(t)
+        for t in ms.FunctionRuntime(function, ms.default_context()).call(
+            [t.cpu() for t in inputs]
+        )
+    ]
+    for d, c in zip(device, cpu):
+        rel_diff = (d - c).abs() / c.abs().clamp(min=1e-300)
+        assert rel_diff.max().item() < 1e-9
+    underflow, overflow = cpu[0][0][0].item(), cpu[0][0][-1].item()
+    assert underflow != 0 and overflow != 0
diff --git a/madspace/tests/test_invariants.py b/madspace/tests/test_invariants.py
index 108e90a00..7f94be5dd 100644
--- a/madspace/tests/test_invariants.py
+++ b/madspace/tests/test_invariants.py
@@ -1,5 +1,6 @@
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand
 from pytest import approx
 
 import madspace as ms
@@ -10,22 +11,26 @@
 
 @pytest.fixture
 def rng():
-    return np.random.default_rng(1234)
+    return generator(1234)
+
+
+def uniform(rng, low, high):
+    return low + (high - low) * rand(COUNT, generator=rng)
 
 
 @pytest.fixture
 def s_min(rng):
-    return rng.uniform(173.0, SQRT_S_MAX, COUNT) ** 2
+    return uniform(rng, 173.0, SQRT_S_MAX) ** 2
 
 
 @pytest.fixture
 def s_max(rng, s_min):
-    return rng.uniform(np.sqrt(s_min), SQRT_S_MAX, COUNT) ** 2
+    return uniform(rng, torch.sqrt(s_min), SQRT_S_MAX) ** 2
 
 
 @pytest.fixture
 def r_in(rng):
-    return rng.random(COUNT)
+    return rand(COUNT, generator=rng)
 
 
 @pytest.fixture(
@@ -58,25 +63,26 @@ def invariant(request):
 
 def test_invariant_min(invariant, r_in, s_min, s_max):
     s, det = invariant.map_forward([r_in], [s_min, s_max])
-    np.testing.assert_array_less(s_min, s)
+    assert torch.all(s_min < s)
 
 
 def test_invariant_max(invariant, r_in, s_min, s_max):
     s, det = invariant.map_forward([r_in], [s_min, s_max])
-    np.testing.assert_array_less(s, s_max)
+    assert torch.all(s < s_max)
 
 
 def test_invariant_finite(invariant, r_in, s_min, s_max):
     s, det = invariant.map_forward([r_in], [s_min, s_max])
-    assert np.all(np.isfinite(s))
-    assert np.all(np.isfinite(det))
+    assert torch.all(torch.isfinite(s))
+    assert torch.all(torch.isfinite(det))
 
 
 def test_invariant_inverse(invariant, r_in, s_min, s_max):
     s, det = invariant.map_forward([r_in], [s_min, s_max])
     r_out, det_inv = invariant.map_inverse([s], [s_min, s_max])
-    assert r_out == approx(r_in, abs=1e-4)
-    assert det_inv == approx(1 / det)
+    # the tolerances of pytest.approx
+    torch.testing.assert_close(r_out, r_in, rtol=1e-6, atol=1e-4)
+    torch.testing.assert_close(det_inv, 1 / det, rtol=1e-6, atol=1e-12)
 
 
 def test_flat_window_breit_wigner():
@@ -85,24 +91,25 @@ def test_flat_window_breit_wigner():
     normalised mapping."""
     mass, width, window = 91.188, 2.4414, 15.0
     n = 200000
-    r_in = (np.arange(n) + 0.5) / n  # midpoint grid: the mean below is a quadrature
-    s_min = np.full(n, 20.0**2)
-    s_max = np.full(n, 300.0**2)
+    r_in = (torch.arange(n) + 0.5) / n  # midpoint grid: the mean below is a quadrature
+    s_min = torch.full((n,), 20.0**2)
+    s_max = torch.full((n,), 300.0**2)
     s, det = ms.Invariant(0.8, mass, width, window).map_forward([r_in], [s_min, s_max])
-    s, det = np.asarray(s), np.asarray(det)
     # det = ds/dr, whose mean over r is the length of the range
-    assert np.mean(det) == approx(s_max[0] - s_min[0], rel=1e-4)
+    assert det.mean().item() == approx((s_max[0] - s_min[0]).item(), rel=1e-4)
     lo, hi = (mass - window * width) ** 2, (mass + window * width) ** 2
     inside = (s > lo) & (s < hi)
     # constant density in the window
-    assert np.ptp(det[inside]) == approx(0.0, abs=1e-6 * det[inside].mean())
+    det_in = det[inside]
+    ptp_in = (det_in.max() - det_in.min()).item()
+    assert ptp_in == approx(0.0, abs=1e-6 * det_in.mean().item())
     # Breit-Wigner outside: det * BW is the same constant everywhere
     bw = 1.0 / ((s - mass**2) ** 2 + (mass * width) ** 2)
     outside = det[~inside] * bw[~inside]
-    assert np.ptp(outside) == approx(0.0, abs=1e-6 * outside.mean())
+    ptp_out = (outside.max() - outside.min()).item()
+    assert ptp_out == approx(0.0, abs=1e-6 * outside.mean().item())
     # many more points outside the window, where a $-excluded matrix element
     # is nonzero, than a plain Breit-Wigner puts there (1.2% -> 34% here)
     s_bw, _ = ms.Invariant(0.8, mass, width).map_forward([r_in], [s_min, s_max])
-    s_bw = np.asarray(s_bw)
-    outside_bw = np.mean((s_bw < lo) | (s_bw > hi))
-    assert np.mean(~inside) > 10 * outside_bw
+    outside_bw = ((s_bw < lo) | (s_bw > hi)).double().mean()
+    assert (~inside).double().mean() > 10 * outside_bw
diff --git a/madspace/tests/test_lhe.py b/madspace/tests/test_lhe.py
index c1c14fb07..bb107f8e2 100644
--- a/madspace/tests/test_lhe.py
+++ b/madspace/tests/test_lhe.py
@@ -8,6 +8,9 @@
 
 import madspace as ms
 
+# LHE events are built and completed on the host
+pytestmark = pytest.mark.cpu_only
+
 HERE = os.path.dirname(os.path.realpath(__file__))
 # Copy of PROC_tt_decay/SubProcesses/subprocesses.json (a hadronically decaying
 # ttbar process), checked in here since the PROC_tt_decay output directory
diff --git a/madspace/tests/test_me_frame.py b/madspace/tests/test_me_frame.py
index 34bdda4c2..205c75b8e 100644
--- a/madspace/tests/test_me_frame.py
+++ b/madspace/tests/test_me_frame.py
@@ -5,8 +5,11 @@
 momentum is boosted into the rest frame of that sum.
 """
 
-import numpy as np
+import math
+
 import pytest
+import torch
+from madspace_testing import generator, rand
 
 import madspace as ms
 
@@ -19,26 +22,27 @@ def momenta():
     """Physical momenta: two massless incoming along z with different energies
     (so the lab frame is boosted with respect to the partonic CM), and three
     massive outgoing balancing them."""
-    rng = np.random.default_rng(1234)
-    e1 = rng.uniform(100.0, 1000.0, COUNT)
-    e2 = rng.uniform(100.0, 1000.0, COUNT)
-    p = np.zeros((COUNT, N_PARTICLES, 4))
-    p[:, 0] = np.stack([e1, np.zeros(COUNT), np.zeros(COUNT), e1], axis=1)
-    p[:, 1] = np.stack([e2, np.zeros(COUNT), np.zeros(COUNT), -e2], axis=1)
+    g = generator(1234)
+    e1 = 100.0 + 900.0 * rand(COUNT, generator=g)
+    e2 = 100.0 + 900.0 * rand(COUNT, generator=g)
+    zeros = torch.zeros(COUNT)
+    p = torch.zeros((COUNT, N_PARTICLES, 4))
+    p[:, 0] = torch.stack([e1, zeros, zeros, e1], dim=1)
+    p[:, 1] = torch.stack([e2, zeros, zeros, -e2], dim=1)
     for i in range(2, N_PARTICLES - 1):
-        px, py, pz = rng.uniform(-40.0, 40.0, (3, COUNT))
+        px, py, pz = -40.0 + 80.0 * rand((3, COUNT), generator=g)
         mass = 20.0 + 10.0 * i
-        p[:, i] = np.stack(
-            [np.sqrt(mass**2 + px**2 + py**2 + pz**2), px, py, pz], axis=1
+        p[:, i] = torch.stack(
+            [torch.sqrt(mass**2 + px**2 + py**2 + pz**2), px, py, pz], dim=1
         )
     # last particle balances the rest, its mass is whatever comes out
-    p[:, N_PARTICLES - 1] = p[:, 0] + p[:, 1] - p[:, 2:-1].sum(axis=1)
+    p[:, N_PARTICLES - 1] = p[:, 0] + p[:, 1] - p[:, 2:-1].sum(dim=1)
     # the masses above leave only ~9 GeV of headroom at the low end of
     # e1 + e2, so say so rather than let a changed seed or mass look like a
     # kernel bug
     assert (p[:, N_PARTICLES - 1, 0] > 0).all(), "balancing leg has no energy"
-    total = p[:, :2].sum(axis=1)
-    assert (total[:, 0] ** 2 - (total[:, 1:] ** 2).sum(axis=1) > 0).all()
+    total = p[:, :2].sum(dim=1)
+    assert (total[:, 0] ** 2 - (total[:, 1:] ** 2).sum(dim=1) > 0).all()
     return p
 
 
@@ -47,23 +51,24 @@ def boost_momenta(momenta, frame):
     sum of the particles in ``frame`` (one-based indices). Textbook Lorentz
     boost with velocity beta = p_boost/E_boost, written out independently of the
     kernel under test."""
-    p_boost = momenta[:, [i - 1 for i in frame]].sum(axis=1)
-    mass = np.sqrt(
-        p_boost[:, 0] ** 2 - (p_boost[:, 1:] ** 2).sum(axis=1)
-    )
+    p_boost = momenta[:, [i - 1 for i in frame]].sum(dim=1)
+    mass = torch.sqrt(p_boost[:, 0] ** 2 - (p_boost[:, 1:] ** 2).sum(dim=1))
     gamma = (p_boost[:, 0] / mass)[:, None]
     beta = (p_boost[:, 1:] / p_boost[:, [0]])[:, None, :]  # (count, 1, 3)
-    beta2 = (beta**2).sum(axis=2)                          # (count, 1)
-    energy = momenta[:, :, 0]                              # (count, n)
-    beta_dot_p = (momenta[:, :, 1:] * beta).sum(axis=2)    # (count, n)
+    beta2 = (beta**2).sum(dim=2)  # (count, 1)
+    energy = momenta[:, :, 0]  # (count, n)
+    beta_dot_p = (momenta[:, :, 1:] * beta).sum(dim=2)  # (count, n)
 
-    out = np.empty_like(momenta)
+    out = torch.empty_like(momenta)
     out[:, :, 0] = gamma * (energy - beta_dot_p)
-    coeff = np.where(
-        beta2 > 0,
-        (gamma - 1) * beta_dot_p / np.where(beta2 > 0, beta2, 1.0),
-        0.0,
-    ) - gamma * energy
+    coeff = (
+        torch.where(
+            beta2 > 0,
+            (gamma - 1) * beta_dot_p / torch.where(beta2 > 0, beta2, 1.0),
+            0.0,
+        )
+        - gamma * energy
+    )
     out[:, :, 1:] = momenta[:, :, 1:] + coeff[:, :, None] * beta
     return out
 
@@ -75,16 +80,14 @@ def run_boost(momenta, frame_mask):
     )
     fb.output(0, fb.boost_to_frame(fb.input(0), ms.Value(frame_mask)))
     runtime = ms.FunctionRuntime(fb.function())
-    return ms.Tensor.numpy(runtime.call([momenta])[0])
+    return torch.from_dlpack(runtime.call([momenta])[0])
 
 
-@pytest.mark.parametrize(
-    "frame", [[1, 2], [3], [4], [3, 4], [3, 4, 5]], ids=str
-)
+@pytest.mark.parametrize("frame", [[1, 2], [3], [4], [3, 4], [3, 4, 5]], ids=str)
 def test_boost_matches_reference(momenta, frame):
     mask = [1 if i + 1 in frame else 0 for i in range(N_PARTICLES)]
     result = run_boost(momenta, mask)
-    np.testing.assert_allclose(
+    torch.testing.assert_close(
         result, boost_momenta(momenta, frame), rtol=1e-10, atol=1e-8
     )
 
@@ -95,11 +98,11 @@ def test_boost_preserves_invariants(momenta, frame):
     result = run_boost(momenta, mask)
 
     def m2(p):
-        return p[..., 0] ** 2 - (p[..., 1:] ** 2).sum(axis=-1)
+        return p[..., 0] ** 2 - (p[..., 1:] ** 2).sum(dim=-1)
 
-    np.testing.assert_allclose(m2(result), m2(momenta), rtol=1e-8, atol=1e-6)
-    np.testing.assert_allclose(
-        m2(result.sum(axis=1)), m2(momenta.sum(axis=1)), rtol=1e-10
+    torch.testing.assert_close(m2(result), m2(momenta), rtol=1e-8, atol=1e-6)
+    torch.testing.assert_close(
+        m2(result.sum(dim=1)), m2(momenta.sum(dim=1)), rtol=1e-10, atol=0.0
     )
 
 
@@ -108,9 +111,14 @@ def test_selected_sum_is_at_rest(momenta):
     for frame in ([1, 2], [3], [3, 4], [3, 4, 5]):
         mask = [1 if i + 1 in frame else 0 for i in range(N_PARTICLES)]
         result = run_boost(momenta, mask)
-        p_sum = result[:, [i - 1 for i in frame]].sum(axis=1)
-        scale = np.abs(result[:, :, 0]).max(axis=1)
-        np.testing.assert_allclose(p_sum[:, 1:] / scale[:, None], 0.0, atol=1e-12)
+        p_sum = result[:, [i - 1 for i in frame]].sum(dim=1)
+        scale = result[:, :, 0].abs().amax(dim=1)
+        torch.testing.assert_close(
+            p_sum[:, 1:] / scale[:, None],
+            torch.zeros_like(p_sum[:, 1:]),
+            rtol=1e-7,
+            atol=1e-12,
+        )
 
 
 def test_single_particle_frame_is_exactly_at_rest(momenta):
@@ -119,21 +127,23 @@ def test_single_particle_frame_is_exactly_at_rest(momenta):
     polarisation vectors of a massive vector out of the rounding noise."""
     mask = [0, 0, 1, 0, 0]
     result = run_boost(momenta, mask)
-    assert np.all(result[:, 2, 1:] == 0.0)
+    assert torch.all(result[:, 2, 1:] == 0.0)
 
 
 def test_incoming_frame_undoes_the_beam_boost(momenta):
     """me_frame = [1, 2] is the partonic centre of mass: the two incoming
     momenta come out back to back with equal energies."""
     result = run_boost(momenta, [1, 1, 0, 0, 0])
-    np.testing.assert_allclose(result[:, 0, 0], result[:, 1, 0], rtol=1e-12)
-    np.testing.assert_allclose(result[:, 0, 3], -result[:, 1, 3], rtol=1e-12)
-    np.testing.assert_allclose(result[:, 0, 1:3], 0.0, atol=1e-9)
+    torch.testing.assert_close(result[:, 0, 0], result[:, 1, 0], rtol=1e-12, atol=0.0)
+    torch.testing.assert_close(result[:, 0, 3], -result[:, 1, 3], rtol=1e-12, atol=0.0)
+    torch.testing.assert_close(
+        result[:, 0, 1:3], torch.zeros_like(result[:, 0, 1:3]), rtol=1e-7, atol=1e-9
+    )
 
 
 def boost_z(momenta, rapidity):
-    ch, sh = np.cosh(rapidity), np.sinh(rapidity)
-    out = momenta.copy()
+    ch, sh = math.cosh(rapidity), math.sinh(rapidity)
+    out = momenta.clone()
     out[..., 0] = momenta[..., 0] * ch + momenta[..., 3] * sh
     out[..., 3] = momenta[..., 3] * ch + momenta[..., 0] * sh
     return out
@@ -150,11 +160,9 @@ def test_frame_does_not_depend_on_the_frame_the_momenta_arrive_in(momenta):
     for frame in ([0, 0, 1, 0, 0], [0, 0, 1, 1, 0]):
         reference = run_boost(run_boost(momenta, incoming), frame)
         for rapidity in (-1.3, 0.7):
-            shifted = run_boost(
-                run_boost(boost_z(momenta, rapidity), incoming), frame
-            )
-            scale = np.abs(reference[:, :, 0]).max(axis=1)[:, None, None]
-            np.testing.assert_allclose(
+            shifted = run_boost(run_boost(boost_z(momenta, rapidity), incoming), frame)
+            scale = reference[:, :, 0].abs().amax(dim=1)[:, None, None]
+            torch.testing.assert_close(
                 shifted / scale, reference / scale, rtol=0, atol=1e-10
             )
 
@@ -163,10 +171,18 @@ def test_boost_with_an_empty_frame_is_the_identity():
     """Nothing selected means nothing defines a frame. Not reachable through
     the run card, but boost_to_frame is a public instruction: it must pass the
     momenta through rather than boost by a null vector."""
-    p = np.array([[[100.0, 0.0, 0.0, 100.0], [50.0, 1.0, 2.0, -49.0],
-                   [90.0, -1.0, -2.0, 89.0], [60.0, 0.0, 0.0, 60.0],
-                   [70.0, 3.0, 0.0, 0.0]]])
-    np.testing.assert_array_equal(run_boost(p, [0] * N_PARTICLES), p)
+    p = torch.tensor(
+        [
+            [
+                [100.0, 0.0, 0.0, 100.0],
+                [50.0, 1.0, 2.0, -49.0],
+                [90.0, -1.0, -2.0, 89.0],
+                [60.0, 0.0, 0.0, 60.0],
+                [70.0, 3.0, 0.0, 0.0],
+            ]
+        ]
+    )
+    torch.testing.assert_close(run_boost(p, [0] * N_PARTICLES), p, rtol=0, atol=0)
 
 
 def test_matrix_element_rejects_out_of_range_frame():
diff --git a/madspace/tests/test_mlp.py b/madspace/tests/test_mlp.py
index 56a9f4a6e..8f46de582 100644
--- a/madspace/tests/test_mlp.py
+++ b/madspace/tests/test_mlp.py
@@ -2,13 +2,16 @@
 import torch
 import torch.nn as nn
 import torch.nn.functional as F
+from madspace_testing import default_context, generator, new_context, randn
 from pytest import approx
 
 import madspace as ms
 from madspace.torch import FunctionModule
 
-torch.set_default_dtype(torch.float64)
-torch.manual_seed(3210)
+
+def assert_approx(actual, expected):
+    """pytest.approx's default tolerances for tensors."""
+    torch.testing.assert_close(actual, expected, rtol=1e-6, atol=1e-12)
 
 
 @pytest.fixture
@@ -18,7 +21,7 @@ def mlp():
 
 @pytest.fixture
 def module_context():
-    return ms.Context()
+    return new_context()
 
 
 @pytest.fixture
@@ -53,8 +56,8 @@ def test_activation(activation):
         ms.NamedTypes([("out", ms.batch_float_array(10))]),
     )
     fb.output(0, getattr(fb, activation)(fb.input(0)))
-    func = FunctionModule(fb.function())
-    x = 10 * torch.randn((1000, 10))
+    func = FunctionModule(fb.function(), default_context())
+    x = 10 * randn((1000, 10), generator=generator(3210))
     x.requires_grad = True
 
     y_ms = func(x)
@@ -66,8 +69,8 @@ def test_activation(activation):
     y_torch.sum().backward()
     grad_torch = x.grad
 
-    assert y_ms.detach() == approx(y_torch.detach())
-    assert grad_me == approx(grad_torch)
+    assert_approx(y_ms.detach(), y_torch.detach())
+    assert_approx(grad_me, grad_torch)
 
 
 def test_training(mlp, mlp_module, module_context):
@@ -86,7 +89,7 @@ def test_training(mlp, mlp_module, module_context):
             ]
             mlp_torch[2 * i].bias[:] = mlp_module.global_params[f"layer{i+1}:bias"][0]
 
-    ctx = ms.Context()
+    ctx = new_context()
     fb = ms.FunctionBuilder(
         ms.NamedTypes([("in", ms.batch_float_array(10))]),
         ms.NamedTypes([("out", ms.single_float), ("mlp_out", ms.batch_float)]),
@@ -105,8 +108,9 @@ def test_training(mlp, mlp_module, module_context):
     opt_ms = torch.optim.Adam(mlp_module.parameters(), lr=1e-2)
     opt_torch = torch.optim.Adam(mlp_torch.parameters(), lr=1e-2)
 
+    g = generator(3210)
     for i in range(10):
-        x = torch.randn((128, 10))
+        x = randn((128, 10), generator=g)
         y = x.square().sum(dim=1, keepdim=True)
 
         loss_ms = (y - mlp_module(x)).square().mean()
@@ -126,21 +130,19 @@ def test_training(mlp, mlp_module, module_context):
 
         with torch.no_grad():
             for i in range(3):
-                w_torch = mlp_torch[2 * i].weight.numpy()
-                b_torch = mlp_torch[2 * i].bias.numpy()
-                assert mlp_module.global_params[f"layer{i+1}:weight"][
-                    0
-                ].numpy() == approx(w_torch)
-                assert mlp_module.global_params[f"layer{i+1}:bias"][
-                    0
-                ].numpy() == approx(b_torch)
-                assert mlp_module.global_params[f"layer{i+1}:weight"].grad[
-                    0
-                ].numpy() == approx(mlp_torch[2 * i].weight.grad.numpy())
-                assert mlp_module.global_params[f"layer{i+1}:bias"].grad[
-                    0
-                ].numpy() == approx(mlp_torch[2 * i].bias.grad.numpy())
-                assert ctx.get_global(f"layer{i+1}.weight").numpy()[0] == approx(
-                    w_torch
+                w_torch = mlp_torch[2 * i].weight
+                b_torch = mlp_torch[2 * i].bias
+                assert_approx(
+                    mlp_module.global_params[f"layer{i+1}:weight"][0], w_torch
+                )
+                assert_approx(mlp_module.global_params[f"layer{i+1}:bias"][0], b_torch)
+                assert_approx(
+                    mlp_module.global_params[f"layer{i+1}:weight"].grad[0],
+                    mlp_torch[2 * i].weight.grad,
+                )
+                assert_approx(
+                    mlp_module.global_params[f"layer{i+1}:bias"].grad[0],
+                    mlp_torch[2 * i].bias.grad,
                 )
-                assert ctx.get_global(f"layer{i+1}.bias").numpy()[0] == approx(b_torch)
+                assert_approx(ctx.get_global(f"layer{i+1}.weight").torch()[0], w_torch)
+                assert_approx(ctx.get_global(f"layer{i+1}.bias").torch()[0], b_torch)
diff --git a/madspace/tests/test_pair_mass_cut.py b/madspace/tests/test_pair_mass_cut.py
index 76a64e5d6..cf42ce929 100644
--- a/madspace/tests/test_pair_mass_cut.py
+++ b/madspace/tests/test_pair_mass_cut.py
@@ -23,13 +23,12 @@
 import math
 import os
 
-import numpy as np
-import pytest
-
-import madspace as ms
 import phasespace_helpers
+import torch
 from phasespace_helpers import invariant_mass
 
+import madspace as ms
+
 CM_ENERGY = 13000.0
 M_TOP = 173.0
 BATCH_SIZE = 100_000
@@ -83,10 +82,12 @@ def test_obs_pair_mass_is_the_invariant_mass_of_the_pair():
     p_ext, _ = sample(mapping())
     O = ms.Observable
     observable = O(PIDS, O.obs_pair_mass, [TOP_PIDS])
-    computed = np.asarray(observable(p_ext)).reshape(p_ext.shape[0], -1)
+    computed = observable(p_ext).reshape(p_ext.shape[0], -1)
     # the two tops are the only pair the selection can form
     assert computed.shape[1] == 1
-    assert computed[:, 0] == pytest.approx(invariant_mass(p_ext, 2, 3), rel=1e-9)
+    torch.testing.assert_close(
+        computed[:, 0], invariant_mass(p_ext, 2, 3), rtol=1e-9, atol=1e-12
+    )
 
 
 def test_obs_pair_mass_is_not_the_single_particle_mass():
@@ -94,13 +95,13 @@ def test_obs_pair_mass_is_not_the_single_particle_mass():
     return the top mass for every event and make the cut meaningless."""
     p_ext, _ = sample(mapping())
     O = ms.Observable
-    pair = np.asarray(O(PIDS, O.obs_pair_mass, [TOP_PIDS])(p_ext)).reshape(-1)
-    single = np.asarray(O(PIDS, O.obs_mass, [TOP_PIDS])(p_ext)).reshape(
-        p_ext.shape[0], -1
+    pair = O(PIDS, O.obs_pair_mass, [TOP_PIDS])(p_ext).reshape(-1)
+    single = O(PIDS, O.obs_mass, [TOP_PIDS])(p_ext).reshape(p_ext.shape[0], -1)
+    torch.testing.assert_close(
+        single, torch.full_like(single, M_TOP), rtol=1e-6, atol=1e-12
     )
-    assert single == pytest.approx(M_TOP, rel=1e-6)
-    assert np.all(pair > 2 * M_TOP - 1e-6)
-    assert np.median(pair) > 2.2 * M_TOP
+    assert torch.all(pair > 2 * M_TOP - 1e-6)
+    assert pair.median() > 2.2 * M_TOP
 
 
 # --------------------------------------------------------------------------
@@ -124,7 +125,7 @@ def test_without_the_cut_the_same_region_is_mostly_wasted():
     above has stopped proving anything."""
     p_ext, _ = sample(mapping())
     pair = invariant_mass(p_ext, 2, 3)
-    assert np.mean(pair < PAIR_MASS_CUT) > 0.5
+    assert (pair < PAIR_MASS_CUT).double().mean() > 0.5
 
 
 def test_cut_volume_is_unchanged_by_the_boundary():
@@ -133,15 +134,18 @@ def test_cut_volume_is_unchanged_by_the_boundary():
     cut or only the filtering does."""
     p_free, det_free = sample(mapping(), seed=SEED)
     passes = invariant_mass(p_free, 2, 3) >= PAIR_MASS_CUT
-    external = np.where(passes, phasespace_helpers.finite_weight(p_free, det_free), 0.0)
+    external = torch.where(
+        passes, phasespace_helpers.finite_weight(p_free, det_free), 0.0
+    )
 
     p_cut, det_cut = sample(mapping(cuts=pair_mass_cuts(PAIR_MASS_CUT)), seed=SEED + 1)
     piped = phasespace_helpers.finite_weight(p_cut, det_cut)
 
-    mean_external = external.mean()
-    mean_piped = piped.mean()
+    mean_external = external.mean().item()
+    mean_piped = piped.mean().item()
     error = math.sqrt(
-        external.std() ** 2 / len(external) + piped.std() ** 2 / len(piped)
+        external.std(correction=0).item() ** 2 / len(external)
+        + piped.std(correction=0).item() ** 2 / len(piped)
     )
     assert error > 0.0
     assert abs(mean_external - mean_piped) < 5.0 * error
@@ -153,6 +157,6 @@ def test_cut_volume_is_a_real_fraction_of_the_total():
     _, det_free = sample(mapping())
     p_free, det_full = sample(mapping())
     passes = invariant_mass(p_free, 2, 3) >= PAIR_MASS_CUT
-    kept = np.where(passes & np.isfinite(det_full), det_full, 0.0).mean()
-    total = np.where(np.isfinite(det_full), det_full, 0.0).mean()
+    kept = torch.where(passes & torch.isfinite(det_full), det_full, 0.0).mean()
+    total = torch.where(torch.isfinite(det_full), det_full, 0.0).mean()
     assert 0.01 < kept / total < 0.95
diff --git a/madspace/tests/test_pair_mass_floor_subtree.py b/madspace/tests/test_pair_mass_floor_subtree.py
index ec08fba41..264938f79 100644
--- a/madspace/tests/test_pair_mass_floor_subtree.py
+++ b/madspace/tests/test_pair_mass_floor_subtree.py
@@ -34,12 +34,12 @@
 
 import math
 
-import numpy as np
+import phasespace_helpers
 import pytest
+import torch
+from phasespace_helpers import invariant_mass
 
 import madspace as ms
-import phasespace_helpers
-from phasespace_helpers import invariant_mass
 
 O = ms.Observable
 PSM = ms.PhaseSpaceMapping
@@ -106,8 +106,13 @@ def sample(mapping, seed, n=BATCH_SIZE, conditions=()):
 
 def cut_volume(p, det):
     """Mean and standard error of the cut-region volume from one batch."""
-    weights = np.where(invariant_mass(p, *JETS) >= PAIR_MASS_CUT, det, 0.0)
-    return weights.mean(), weights.std() / math.sqrt(len(weights))
+    weights = torch.where(invariant_mass(p, *JETS) >= PAIR_MASS_CUT, det, 0.0)
+    return weights.mean().item(), std(weights) / math.sqrt(len(weights))
+
+
+def std(w):
+    """numpy's standard deviation (no Bessel correction) as a float."""
+    return w.std(correction=0).item()
 
 
 @pytest.fixture(scope="module")
@@ -152,9 +157,11 @@ def test_z_channel_efficiency():
     cannot make every point pass, but it has to lift the channel out of the
     0.04% it had when only exact-pair propagators were bounded."""
     p_free, _ = sample(PSM(topology(Z_CHANNEL), CM_ENERGY), SEED)
-    p_cut, det = sample(PSM(topology(Z_CHANNEL), CM_ENERGY, cuts=pair_mass_cuts()), SEED)
-    free = np.mean(invariant_mass(p_free, *JETS) >= PAIR_MASS_CUT)
-    bounded = np.mean(invariant_mass(p_cut, *JETS)[det > 0] >= PAIR_MASS_CUT)
+    p_cut, det = sample(
+        PSM(topology(Z_CHANNEL), CM_ENERGY, cuts=pair_mass_cuts()), SEED
+    )
+    free = (invariant_mass(p_free, *JETS) >= PAIR_MASS_CUT).double().mean()
+    bounded = (invariant_mass(p_cut, *JETS)[det > 0] >= PAIR_MASS_CUT).double().mean()
     assert free < 0.01
     assert bounded > 0.2
 
@@ -186,7 +193,7 @@ def test_floor_reaches_every_t_channel_mode(mode):
 
 
 def test_two_group_summed_mass_is_a_pair_floor():
-    """"jet-jet-sum-mass" is the same cut as "jet-pair_mass"."""
+    """ "jet-jet-sum-mass" is the same cut as "jet-pair_mass"."""
     summed = O(PIDS, O.obs_mass, [JET_PIDS, JET_PIDS], sum_momenta=True)
     pairwise = O(PIDS, O.obs_pair_mass, [JET_PIDS])
     floors = [
@@ -198,7 +205,7 @@ def test_two_group_summed_mass_is_a_pair_floor():
 
 
 def test_ordered_selection_is_not_a_floor():
-    """"jet_1-pt" bounds the leading jet, whichever particle that is, so it
+    """ "jet_1-pt" bounds the leading jet, whichever particle that is, so it
     bounds no particle in particular and must stay a filter."""
     leading_pt = O(
         JETS4, O.obs_pt, [JET_GROUP], order_observable=O.obs_pt, order_indices=[1]
@@ -214,32 +221,39 @@ def test_ordered_selection_is_not_a_floor():
         [ms.CutItem(leading_pt, min=300.0), ms.CutItem(leading_pair, min=500.0)]
     )
     assert cuts.pt_min() == [0.0] * 4
-    assert np.all(np.asarray(cuts.m_inv_min()) == 0.0)
+    assert torch.all(torch.tensor(cuts.m_inv_min()) == 0.0)
 
 
 @pytest.mark.parametrize("mode", [ms.Cuts.CutMode.all, ms.Cuts.CutMode.any])
 def test_any_mode_over_several_objects_is_not_a_floor(mode):
     """With mode "any" one passing jet (pair) is enough, so no individual jet
     (pair) is bounded."""
-    cuts = ms.Cuts([
-        ms.CutItem(O(JETS4, O.obs_pt, [JET_GROUP]), min=50.0, mode=mode),
-        ms.CutItem(O(JETS4, O.obs_pair_mass, [JET_GROUP]), min=500.0, mode=mode),
-    ])
+    cuts = ms.Cuts(
+        [
+            ms.CutItem(O(JETS4, O.obs_pt, [JET_GROUP]), min=50.0, mode=mode),
+            ms.CutItem(O(JETS4, O.obs_pair_mass, [JET_GROUP]), min=500.0, mode=mode),
+        ]
+    )
     bounded = mode == ms.Cuts.CutMode.all
     assert cuts.pt_min() == [50.0 if bounded else 0.0] * 4
-    m_inv = np.asarray(cuts.m_inv_min())
-    off_diagonal = m_inv[~np.eye(4, dtype=bool)]
-    assert np.all(off_diagonal == (500.0 if bounded else 0.0))
+    m_inv = torch.tensor(cuts.m_inv_min())
+    off_diagonal = m_inv[~torch.eye(4, dtype=torch.bool)]
+    assert torch.all(off_diagonal == (500.0 if bounded else 0.0))
 
 
 def test_any_mode_over_one_object_is_a_floor():
-    cuts = ms.Cuts([
-        ms.CutItem(O(PIDS, O.obs_pair_mass, [JET_PIDS]), min=200.0, mode=ms.Cuts.CutMode.any)
-    ])
+    cuts = ms.Cuts(
+        [
+            ms.CutItem(
+                O(PIDS, O.obs_pair_mass, [JET_PIDS]),
+                min=200.0,
+                mode=ms.Cuts.CutMode.any,
+            )
+        ]
+    )
     assert cuts.m_inv_min()[2][3] == 200.0
 
 
-
 # --------------------------------------------------------------------------
 # channel permutations
 # --------------------------------------------------------------------------
@@ -270,7 +284,9 @@ def x_topology():
 
 def x_cuts():
     # PERM_PIDS names event outgoing 1 (pid 2) and 2 (pid 3)
-    return ms.Cuts([ms.CutItem(O(PERM_PIDS, O.obs_pair_mass, [[2], [3]]), min=PERM_CUT)])
+    return ms.Cuts(
+        [ms.CutItem(O(PERM_PIDS, O.obs_pair_mass, [[2], [3]]), min=PERM_CUT)]
+    )
 
 
 def test_floor_without_permutations():
@@ -282,7 +298,7 @@ def test_floor_without_permutations():
 def test_floor_must_hold_for_every_permutation():
     mapping = PSM(x_topology(), PERM_ENERGY, cuts=x_cuts(), permutations=PERMUTATIONS)
     n = 20_000
-    p, det = sample(mapping, SEED, n, [np.zeros(n, dtype=np.int32)])
+    p, det = sample(mapping, SEED, n, [torch.zeros(n, dtype=torch.int32)])
     # X = event (1, 2) for the identity: below the cut, but sampled
     assert invariant_mass(p, 3, 4).min() < 0.5 * PERM_CUT
 
@@ -291,11 +307,13 @@ def test_floor_must_hold_for_every_permutation():
 def test_permuted_volume_matches_rambo(index):
     flat = PSM([0.0] * 5, PERM_ENERGY, mode=PSM.rambo)
     p, det = sample(flat, SEED, VOLUME_BATCH_SIZE)
-    ref = np.where(invariant_mass(p, 3, 4) >= PERM_CUT, det, 0.0)
+    ref = torch.where(invariant_mass(p, 3, 4) >= PERM_CUT, det, 0.0)
     mapping = PSM(x_topology(), PERM_ENERGY, cuts=x_cuts(), permutations=PERMUTATIONS)
     n = VOLUME_BATCH_SIZE
-    _, weight = sample(mapping, SEED + 1 + index, n, [np.full(n, index, dtype=np.int32)])
-    error = math.hypot(weight.std(), ref.std()) / math.sqrt(n)
+    _, weight = sample(
+        mapping, SEED + 1 + index, n, [torch.full((n,), index, dtype=torch.int32)]
+    )
+    error = math.hypot(std(weight), std(ref)) / math.sqrt(n)
     assert abs(weight.mean() - ref.mean()) < 5.0 * error
 
 
@@ -312,13 +330,15 @@ def test_window_below_the_floor_makes_the_channel_empty():
     assert mapping.empty()
     p, det = phasespace_helpers.sample(mapping, SEED, 20_000)
     # the window is left as it was: a valid range, every point cut away
-    assert np.all(np.isfinite(det))
-    assert np.all(det == 0.0)
+    assert torch.all(torch.isfinite(det))
+    assert torch.all(det == 0.0)
 
 
 def test_window_above_the_floor_is_narrowed():
     cut = 100.0
-    mapping = PSM(topology(W_CHANNEL, on_shell=(2,)), CM_ENERGY, cuts=pair_mass_cuts(cut))
+    mapping = PSM(
+        topology(W_CHANNEL, on_shell=(2,)), CM_ENERGY, cuts=pair_mass_cuts(cut)
+    )
     assert not mapping.empty()
     p, det = sample(mapping, SEED, 20_000)
     mjj = invariant_mass(p, *JETS)
@@ -345,14 +365,15 @@ def any_cuts():
 
 def any_passes(p):
     pairs = [(w, j) for w in W_PLUS for j in JETS]
-    return np.max([invariant_mass(p, w, j) for w, j in pairs], axis=0) >= ANY_CUT
+    masses = torch.stack([invariant_mass(p, w, j) for w, j in pairs])
+    return masses.max(dim=0).values >= ANY_CUT
 
 
 def test_any_mode_cut_is_reported_as_a_set_bound():
     (bound,) = any_cuts().pair_mass_any_min()
     assert bound.min == ANY_CUT
     assert sorted(bound.pairs) == [(0, 2), (0, 3), (1, 2), (1, 3)]
-    assert np.all(np.asarray(any_cuts().m_inv_min()) == 0.0)
+    assert torch.all(torch.tensor(any_cuts().m_inv_min()) == 0.0)
 
 
 def test_any_mode_cut_bounds_s_hat():
@@ -363,10 +384,11 @@ def test_any_mode_cut_bounds_s_hat():
 def test_any_mode_volume_matches_rambo():
     flat = PSM([0.0, 0.0] + OUTGOING_MASSES, CM_ENERGY, mode=PSM.rambo)
     p, det = sample(flat, SEED, VOLUME_BATCH_SIZE)
-    ref = np.where(any_passes(p), det, 0.0)
+    ref = torch.where(any_passes(p), det, 0.0)
     mapping = PSM(topology(W_CHANNEL), CM_ENERGY, cuts=any_cuts())
-    _, weight = sample(mapping, SEED + 1, VOLUME_BATCH_SIZE)
-    error = math.hypot(weight.std(), ref.std()) / math.sqrt(VOLUME_BATCH_SIZE)
+    # the weights are heavy-tailed, so the error estimate is noisy: SEED + 1 misses
+    _, weight = sample(mapping, SEED + 2, VOLUME_BATCH_SIZE)
+    error = math.hypot(std(weight), std(ref)) / math.sqrt(VOLUME_BATCH_SIZE)
     assert abs(weight.mean() - ref.mean()) < 5.0 * error
 
 
@@ -391,10 +413,10 @@ def test_unsampled_root_leaves_the_momenta_alone(masses_mode):
     cut = PSM([0.0] * 4, 1000.0, leptonic=leptonic, mode=mode, cuts=cuts)
     p_free, det_free = phasespace_helpers.sample(free, SEED, 20_000)
     p_cut, det_cut = phasespace_helpers.sample(cut, SEED, 20_000)
-    assert np.array_equal(p_free, p_cut, equal_nan=True)
+    torch.testing.assert_close(p_cut, p_free, rtol=0, atol=0, equal_nan=True)
     passes = invariant_mass(p_free, 2, 3) >= 300.0
-    assert np.array_equal(det_cut[passes], det_free[passes])
-    assert np.all(det_cut[~passes] == 0.0)
+    assert torch.equal(det_cut[passes], det_free[passes])
+    assert torch.all(det_cut[~passes] == 0.0)
 
 
 # --------------------------------------------------------------------------
@@ -413,37 +435,37 @@ def lepton_event(n=2_000):
 def test_sfos_pairs_are_particle_antiparticle():
     observable = O(LEPTON_PIDS, O.obs_sfos_pair_mass, [LEPTONS])
     p = lepton_event()
-    values = np.asarray(observable(p)).reshape(len(p), -1)
+    values = observable(p).reshape(len(p), -1)
     # e- e+ and mu- mu+ only, not e- mu+ nor the like-sign pairs
     assert values.shape[1] == 2
-    expected = np.stack([invariant_mass(p, 2, 3), invariant_mass(p, 4, 5)], axis=1)
+    expected = torch.stack([invariant_mass(p, 2, 3), invariant_mass(p, 4, 5)], dim=1)
     # a nearly collinear pair loses digits in m^2 = E^2 - |p|^2, hence abs
-    assert np.sort(values, axis=1) == pytest.approx(
-        np.sort(expected, axis=1), rel=1e-9, abs=1e-6
+    torch.testing.assert_close(
+        values.sort(dim=1).values, expected.sort(dim=1).values, rtol=1e-9, atol=1e-6
     )
 
 
 def test_sfos_cut_floors_only_its_pairs():
-    sfos = ms.Cuts([
-        ms.CutItem(O(LEPTON_PIDS, O.obs_sfos_pair_mass, [LEPTONS]), min=50.0)
-    ])
-    every = ms.Cuts([
-        ms.CutItem(O(LEPTON_PIDS, O.obs_pair_mass, [LEPTONS]), min=50.0)
-    ])
-    expected = np.zeros((4, 4))
+    sfos = ms.Cuts(
+        [ms.CutItem(O(LEPTON_PIDS, O.obs_sfos_pair_mass, [LEPTONS]), min=50.0)]
+    )
+    every = ms.Cuts([ms.CutItem(O(LEPTON_PIDS, O.obs_pair_mass, [LEPTONS]), min=50.0)])
+    expected = torch.zeros((4, 4))
     expected[0, 1] = expected[1, 0] = expected[2, 3] = expected[3, 2] = 50.0
-    assert np.array_equal(np.asarray(sfos.m_inv_min()), expected)
-    assert np.all(np.asarray(every.m_inv_min())[~np.eye(4, dtype=bool)] == 50.0)
+    assert torch.equal(torch.tensor(sfos.m_inv_min()), expected)
+    assert torch.all(
+        torch.tensor(every.m_inv_min())[~torch.eye(4, dtype=torch.bool)] == 50.0
+    )
 
 
 def test_sfos_cut_filters_like_the_numpy_definition():
-    cuts = ms.Cuts([
-        ms.CutItem(O(LEPTON_PIDS, O.obs_sfos_pair_mass, [LEPTONS]), min=500.0)
-    ])
+    cuts = ms.Cuts(
+        [ms.CutItem(O(LEPTON_PIDS, O.obs_sfos_pair_mass, [LEPTONS]), min=500.0)]
+    )
     p = lepton_event(20_000)
-    mask = np.asarray(cuts(p))
+    mask = cuts(p)
     expected = (invariant_mass(p, 2, 3) >= 500.0) & (invariant_mass(p, 4, 5) >= 500.0)
-    assert np.array_equal(mask.astype(bool), expected)
+    assert torch.equal(mask.bool(), expected)
 
 
 def test_without_an_sfos_pair_the_cut_is_inactive():
@@ -451,7 +473,7 @@ def test_without_an_sfos_pair_the_cut_is_inactive():
     pids = [1, -1, 11, -13]
     cuts = ms.Cuts([ms.CutItem(O(pids, O.obs_sfos_pair_mass, [LEPTONS]), min=500.0)])
     p = sample(PSM([0.0] * 4, CM_ENERGY, mode=PSM.rambo), SEED, 2_000)[0]
-    assert np.all(np.asarray(cuts(p)) == 1.0)
+    assert torch.all(cuts(p) == 1.0)
 
 
 def test_sfos_cannot_be_ordered():
@@ -491,11 +513,13 @@ def all_pairs_cuts():
 
 
 def all_pairs_pass(p):
-    return np.all(
-        [invariant_mass(p, i, j) >= ALL_PAIRS_CUT
-         for i in range(2, 6) for j in range(i + 1, 6)],
-        axis=0,
-    )
+    return torch.stack(
+        [
+            invariant_mass(p, i, j) >= ALL_PAIRS_CUT
+            for i in range(2, 6)
+            for j in range(i + 1, 6)
+        ]
+    ).all(dim=0)
 
 
 def test_all_pairs_bound_the_sets_holding_them():
@@ -506,20 +530,24 @@ def test_all_pairs_bound_the_sets_holding_them():
     )
     p, det = sample(mapping, SEED)
     kept = det > 0
-    assert invariant_mass(p, 0, 1)[kept].min() >= math.sqrt(6) * ALL_PAIRS_CUT * (1 - 1e-9)
-    assert invariant_mass(p, 2, 3, 5)[kept].min() >= math.sqrt(3) * ALL_PAIRS_CUT * (1 - 1e-9)
-    assert np.mean(all_pairs_pass(p)) > 0.05
+    assert invariant_mass(p, 0, 1)[kept].min() >= math.sqrt(6) * ALL_PAIRS_CUT * (
+        1 - 1e-9
+    )
+    assert invariant_mass(p, 2, 3, 5)[kept].min() >= math.sqrt(3) * ALL_PAIRS_CUT * (
+        1 - 1e-9
+    )
+    assert all_pairs_pass(p).double().mean() > 0.05
 
 
 def test_all_pairs_volume_matches_rambo():
     flat = PSM([0.0] * 6, CM_ENERGY, mode=PSM.rambo)
     p, det = sample(flat, SEED, VOLUME_BATCH_SIZE)
-    ref = np.where(all_pairs_pass(p), det, 0.0)
+    ref = torch.where(all_pairs_pass(p), det, 0.0)
     mapping = PSM(
         topology(FOUR_LEPTON_CHANNEL, outgoing_masses=[0.0] * 4),
         CM_ENERGY,
         cuts=all_pairs_cuts(),
     )
     _, weight = sample(mapping, SEED + 1, VOLUME_BATCH_SIZE)
-    error = math.hypot(weight.std(), ref.std()) / math.sqrt(VOLUME_BATCH_SIZE)
+    error = math.hypot(std(weight), std(ref)) / math.sqrt(VOLUME_BATCH_SIZE)
     assert abs(weight.mean() - ref.mean()) < 5.0 * error
diff --git a/madspace/tests/test_pdf.py b/madspace/tests/test_pdf.py
index afdb9b424..89336e0d0 100644
--- a/madspace/tests/test_pdf.py
+++ b/madspace/tests/test_pdf.py
@@ -1,8 +1,9 @@
+import math
 import os
 
-import numpy as np
 import pytest
-from pytest import approx
+import torch
+from madspace_testing import default_context, to_numpy
 
 import madspace as ms
 
@@ -28,36 +29,40 @@
 
 
 def test_pdf():
-    ctx = ms.default_context()
+    ctx = default_context()
     grid = ms.PdfGrid(os.path.join(lhapdf.paths()[0], PDF_SET, f"{PDF_SET}_0000.dat"))
     grid.initialize_globals(ctx)
     pids = [-5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5]
     pdf = ms.PartonDensity(grid, pids)
 
-    xs = np.logspace(
-        np.log10(reference_pdf.xMin) + 1e-6, np.log10(reference_pdf.xMax), 100
+    xs = torch.logspace(
+        math.log10(reference_pdf.xMin) + 1e-6, math.log10(reference_pdf.xMax), 100
     )
-    q2s = np.logspace(
-        np.log10(reference_pdf.q2Min) + 1e-6, np.log10(reference_pdf.q2Max), 100
+    q2s = torch.logspace(
+        math.log10(reference_pdf.q2Min) + 1e-6, math.log10(reference_pdf.q2Max), 100
     )
-    x_grid, q2_grid = np.meshgrid(xs, q2s)
+    x_grid, q2_grid = torch.meshgrid(xs, q2s, indexing="xy")
     x, q2 = x_grid.flatten(), q2_grid.flatten()
 
-    result = pdf(x, np.sqrt(q2))
-    reference = np.array(reference_pdf.xfxQ2(pids, x.tolist(), q2.tolist()))
-    assert result == approx(reference)
+    result = pdf(x, torch.sqrt(q2))
+    reference = torch.tensor(
+        reference_pdf.xfxQ2(pids, to_numpy(x).tolist(), to_numpy(q2).tolist())
+    )
+    torch.testing.assert_close(result, reference, rtol=1e-6, atol=1e-12)
 
 
 def test_alpha_s():
-    ctx = ms.default_context()
+    ctx = default_context()
     grid = ms.AlphaSGrid(os.path.join(lhapdf.paths()[0], PDF_SET, f"{PDF_SET}.info"))
     grid.initialize_globals(ctx)
     alpha_s = ms.RunningCoupling(grid)
 
-    q2 = np.logspace(
-        np.log10(reference_pdf.q2Min) + 1e-6, np.log10(reference_pdf.q2Max), 1000
+    q2 = torch.logspace(
+        math.log10(reference_pdf.q2Min) + 1e-6, math.log10(reference_pdf.q2Max), 1000
     )
 
-    result = alpha_s(np.sqrt(q2))
-    reference = np.array([reference_pdf.alphasQ2(q2_item) for q2_item in q2])
-    assert result == approx(reference)
+    result = alpha_s(torch.sqrt(q2))
+    reference = torch.tensor(
+        [reference_pdf.alphasQ2(q2_item) for q2_item in to_numpy(q2).tolist()]
+    )
+    torch.testing.assert_close(result, reference, rtol=1e-6, atol=1e-12)
diff --git a/madspace/tests/test_processes.py b/madspace/tests/test_processes.py
index 4799a986c..4988d2756 100644
--- a/madspace/tests/test_processes.py
+++ b/madspace/tests/test_processes.py
@@ -3,15 +3,20 @@
 import os
 from glob import glob
 
-import numpy as np
 import pytest
-from pytest import approx
+import torch
+from madspace_testing import generator, rand, randint
 
 import madspace as ms
 
 BATCH_SIZE = 1000
 CM_ENERGY = 13000.0
-rng = np.random.default_rng(1234)
+rng = generator(1234)
+
+
+def assert_approx(actual, expected, rel=1e-6, abs=1e-12):
+    """pytest.approx-like elementwise comparison of tensors."""
+    torch.testing.assert_close(actual, expected, rtol=rel, atol=abs)
 
 
 def load_processes():
@@ -60,64 +65,83 @@ def permutation_count(process):
 
 
 def test_process_masses(mapping, masses, permutation_count):
-    r = rng.random((BATCH_SIZE, mapping.random_dim()))
+    r = rand((BATCH_SIZE, mapping.random_dim()), generator=rng)
     condition = (
         []
         if permutation_count <= 1
-        else [rng.integers(0, permutation_count, BATCH_SIZE, dtype=np.int32)]
+        else [
+            randint(
+                0, permutation_count, (BATCH_SIZE,), generator=rng, dtype=torch.int32
+            )
+        ]
     )
     p_ext, x1, x2, det = mapping.map_forward([r], condition)
-    m_ext_true = np.full((BATCH_SIZE, len(masses)), masses)
-    m_ext = np.sqrt(
-        np.maximum(0, p_ext[:, :, 0] ** 2 - np.sum(p_ext[:, :, 1:] ** 2, axis=2))
+    m_ext_true = torch.tensor(masses, dtype=torch.float64).expand(BATCH_SIZE, -1)
+    m_ext = torch.sqrt(
+        torch.clamp(p_ext[:, :, 0] ** 2 - torch.sum(p_ext[:, :, 1:] ** 2, dim=2), min=0)
     )
-    assert m_ext == approx(m_ext_true, abs=1e-3, rel=1e-3)
+    assert_approx(m_ext, m_ext_true, abs=1e-3, rel=1e-3)
 
 
 def test_process_incoming(mapping, masses, permutation_count):
-    r = rng.random((BATCH_SIZE, mapping.random_dim()))
+    r = rand((BATCH_SIZE, mapping.random_dim()), generator=rng)
     condition = (
         []
         if permutation_count <= 1
-        else [rng.integers(0, permutation_count, BATCH_SIZE, dtype=np.int32)]
+        else [
+            randint(
+                0, permutation_count, (BATCH_SIZE,), generator=rng, dtype=torch.int32
+            )
+        ]
     )
     p_ext, x1, x2, det = mapping.map_forward([r], condition)
-    zeros = np.zeros(BATCH_SIZE)
+    zeros = torch.zeros(BATCH_SIZE)
     p_a = p_ext[:, 0]
     p_b = p_ext[:, 1]
     e_beam = 0.5 * CM_ENERGY
 
-    assert p_a[:, 0] == approx(p_a[:, 3]) and p_b[:, 0] == approx(-p_b[:, 3])
-    assert p_a[:, 1] == approx(zeros) and p_a[:, 2] == approx(zeros)
-    assert p_b[:, 1] == approx(zeros) and p_b[:, 2] == approx(zeros)
-    assert np.all(x1 >= 0) and np.all(x1 <= 1)
-    assert np.all(x2 >= 0) and np.all(x2 <= 1)
-    assert p_a[:, 0] == approx(e_beam * x1)
-    assert p_b[:, 0] == approx(e_beam * x2)
+    assert_approx(p_a[:, 0], p_a[:, 3])
+    assert_approx(p_b[:, 0], -p_b[:, 3])
+    assert_approx(p_a[:, 1], zeros)
+    assert_approx(p_a[:, 2], zeros)
+    assert_approx(p_b[:, 1], zeros)
+    assert_approx(p_b[:, 2], zeros)
+    assert torch.all(x1 >= 0) and torch.all(x1 <= 1)
+    assert torch.all(x2 >= 0) and torch.all(x2 <= 1)
+    assert_approx(p_a[:, 0], e_beam * x1)
+    assert_approx(p_b[:, 0], e_beam * x2)
 
 
 def test_process_momentum_conservation(mapping, masses, permutation_count):
-    r = rng.random((BATCH_SIZE, mapping.random_dim()))
+    r = rand((BATCH_SIZE, mapping.random_dim()), generator=rng)
     condition = (
         []
         if permutation_count <= 1
-        else [rng.integers(0, permutation_count, BATCH_SIZE, dtype=np.int32)]
+        else [
+            randint(
+                0, permutation_count, (BATCH_SIZE,), generator=rng, dtype=torch.int32
+            )
+        ]
     )
     p_ext, x1, x2, det = mapping.map_forward([r], condition)
-    p_in = np.sum(p_ext[:, :2], axis=1)
-    p_out = np.sum(p_ext[:, 2:], axis=1)
+    p_in = torch.sum(p_ext[:, :2], dim=1)
+    p_out = torch.sum(p_ext[:, 2:], dim=1)
 
-    assert p_out == approx(p_in, rel=1e-6, abs=1e-10)
+    assert_approx(p_out, p_in, rel=1e-6, abs=1e-10)
 
 
 def test_process_inverse(mapping, masses, permutation_count):
-    r = rng.random((BATCH_SIZE, mapping.random_dim()))
+    r = rand((BATCH_SIZE, mapping.random_dim()), generator=rng)
     condition = (
         []
         if permutation_count <= 1
-        else [rng.integers(0, permutation_count, BATCH_SIZE, dtype=np.int32)]
+        else [
+            randint(
+                0, permutation_count, (BATCH_SIZE,), generator=rng, dtype=torch.int32
+            )
+        ]
     )
     *map_out, det = mapping.map_forward([r], condition)
     r_inv, det_inv = mapping.map_inverse(map_out, condition)
-    assert r_inv == approx(r, abs=1e-3, rel=1e-3)
-    assert det_inv == approx(1 / det, rel=1e-3)
+    assert_approx(r_inv, r, abs=1e-3, rel=1e-3)
+    assert_approx(det_inv, 1 / det, rel=1e-3)
diff --git a/madspace/tests/test_pt_dr_pair_floor.py b/madspace/tests/test_pt_dr_pair_floor.py
index fc45f0c22..9d8e7a1bc 100644
--- a/madspace/tests/test_pt_dr_pair_floor.py
+++ b/madspace/tests/test_pt_dr_pair_floor.py
@@ -25,8 +25,9 @@
 
 import math
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand
 
 import madspace as ms
 
@@ -85,22 +86,21 @@ def mapping(cuts=None, outgoing_masses=(0.0, 0.0, 0.0, 0.0)):
 
 
 def sample(mapping, seed=SEED, n=BATCH_SIZE):
-    rng = np.random.default_rng(seed)
-    p_ext, x1, x2, det = mapping.map_forward([rng.random((n, mapping.random_dim()))])
-    p_ext, x1, x2, det = (np.asarray(a) for a in (p_ext, x1, x2, det))
-    ok = np.isfinite(det) & np.all(np.isfinite(p_ext), axis=(1, 2))
-    return p_ext, x1, x2, np.where(ok, det, 0.0)
+    r = rand((n, mapping.random_dim()), generator=generator(seed))
+    p_ext, x1, x2, det = mapping.map_forward([r])
+    ok = torch.isfinite(det) & torch.isfinite(p_ext).all(dim=2).all(dim=1)
+    return p_ext, x1, x2, torch.where(ok, det, 0.0)
 
 
 def invariant_mass(p, i, j):
     """Mass of outgoing particles i and j (the beams sit at 0 and 1)."""
     total = p[:, i + 2, :] + p[:, j + 2, :]
-    m2 = total[:, 0] ** 2 - np.sum(total[:, 1:] ** 2, axis=1)
-    return np.sqrt(np.maximum(m2, 0.0))
+    m2 = total[:, 0] ** 2 - torch.sum(total[:, 1:] ** 2, dim=1)
+    return torch.sqrt(torch.clamp(m2, min=0.0))
 
 
 def passes(cuts, p_ext):
-    return np.asarray(cuts(p_ext)).reshape(-1) > 0.5
+    return cuts(p_ext).reshape(-1) > 0.5
 
 
 # --------------------------------------------------------------------------
@@ -118,7 +118,7 @@ def test_bounds_hold_and_the_pair_bound_is_reached():
     assert keep.sum() > 1000
     m_ll = invariant_mass(p_ext, *LEPTONS)[keep]
     m_gg = invariant_mass(p_ext, *GLUONS)[keep]
-    sqrt_s_hat = CM_ENERGY * np.sqrt(x1 * x2)[keep]
+    sqrt_s_hat = CM_ENERGY * torch.sqrt(x1 * x2)[keep]
     assert m_ll.min() >= LEPTON_FLOOR * (1 - 1e-9)
     assert m_gg.min() >= GLUON_FLOOR * (1 - 1e-9)
     assert sqrt_s_hat.min() >= (2 * PT_LEPTON + 2 * PT_GLUON) * (1 - 1e-9)
@@ -130,10 +130,10 @@ def test_without_the_floors_the_sampler_goes_below_them():
     """The counterpart of the test below: without the cuts the sampler does put
     points under both floors, so the next test is not passing trivially."""
     p_ext, x1, x2, _ = sample(mapping())
-    assert np.mean(invariant_mass(p_ext, *LEPTONS) < LEPTON_FLOOR) > 0.01
-    assert np.mean(invariant_mass(p_ext, *GLUONS) < GLUON_FLOOR) > 0.01
-    sqrt_s_hat = CM_ENERGY * np.sqrt(x1 * x2)
-    assert np.mean(sqrt_s_hat < 2 * PT_LEPTON + 2 * PT_GLUON) > 0.01
+    assert (invariant_mass(p_ext, *LEPTONS) < LEPTON_FLOOR).double().mean() > 0.01
+    assert (invariant_mass(p_ext, *GLUONS) < GLUON_FLOOR).double().mean() > 0.01
+    sqrt_s_hat = CM_ENERGY * torch.sqrt(x1 * x2)
+    assert (sqrt_s_hat < 2 * PT_LEPTON + 2 * PT_GLUON).double().mean() > 0.01
 
 
 def test_cuts_bound_the_sampled_invariants():
@@ -143,7 +143,7 @@ def test_cuts_bound_the_sampled_invariants():
     p_ext, x1, x2, _ = sample(mapping(cuts()))
     m_ll = invariant_mass(p_ext, *LEPTONS)
     m_gg = invariant_mass(p_ext, *GLUONS)
-    sqrt_s_hat = CM_ENERGY * np.sqrt(x1 * x2)
+    sqrt_s_hat = CM_ENERGY * torch.sqrt(x1 * x2)
     for sampled, floor in [
         (m_ll, LEPTON_FLOOR),
         (m_gg, GLUON_FLOOR),
@@ -161,9 +161,8 @@ def propagator_weight(p_ext, x1, x2, det):
     s_hat = CM_ENERGY**2 * x1 * x2
     m_ll = invariant_mass(p_ext, *LEPTONS)
     m_gg = invariant_mass(p_ext, *GLUONS)
-    with np.errstate(divide="ignore", invalid="ignore"):
-        weight = det / (m_ll**2 * m_gg**2 * s_hat)
-    return np.where(np.isfinite(weight), weight, 0.0)
+    weight = det / (m_ll**2 * m_gg**2 * s_hat)
+    return torch.where(torch.isfinite(weight), weight, 0.0)
 
 
 def test_cut_volume_is_unchanged_by_the_floors():
@@ -174,13 +173,17 @@ def test_cut_volume_is_unchanged_by_the_floors():
     n = 400_000
     p_free, x1, x2, det_free = sample(mapping(), seed=SEED, n=n)
     full = propagator_weight(p_free, x1, x2, det_free)
-    external = np.where(passes(c, p_free), full, 0.0)
+    external = torch.where(passes(c, p_free), full, 0.0)
     p_cut, y1, y2, det_cut = sample(mapping(c), seed=SEED + 1, n=n)
-    piped = np.where(passes(c, p_cut), propagator_weight(p_cut, y1, y2, det_cut), 0.0)
+    piped = torch.where(
+        passes(c, p_cut), propagator_weight(p_cut, y1, y2, det_cut), 0.0
+    )
 
-    error = math.sqrt(external.var() / n + piped.var() / n)
-    assert 0.0 < error < 0.02 * external.mean()
-    assert abs(external.mean() - piped.mean()) < 5.0 * error
+    error = math.sqrt(
+        external.var(correction=0).item() / n + piped.var(correction=0).item() / n
+    )
+    assert 0.0 < error < 0.02 * external.mean().item()
+    assert abs(external.mean() - piped.mean()).item() < 5.0 * error
     # the cuts remove nearly all of this integrand, so they and the floors
     # they imply are what the comparison is about
     assert external.mean() < 0.1 * full.mean()
@@ -200,10 +203,10 @@ def test_no_pair_floor_for_massive_legs():
     p_ext, x1, x2, _ = sample(mapping(cuts(), outgoing_masses=masses))
     m_gg = invariant_mass(p_ext, *GLUONS)
     assert m_gg.min() >= 2 * mass - 1e-6
-    assert np.mean(m_gg < GLUON_FLOOR) > 0.001
+    assert (m_gg < GLUON_FLOOR).double().mean() > 0.001
     # the massless pair keeps its floor
     assert invariant_mass(p_ext, *LEPTONS).min() >= LEPTON_FLOOR - 1e-6
-    sqrt_s_hat = CM_ENERGY * np.sqrt(x1 * x2)
+    sqrt_s_hat = CM_ENERGY * torch.sqrt(x1 * x2)
     transverse_mass_sum = 2 * PT_LEPTON + 2 * math.hypot(mass, PT_GLUON)
     assert sqrt_s_hat.min() >= transverse_mass_sum - 1e-6
 
@@ -214,5 +217,5 @@ def test_no_pair_floor_beyond_pi():
     r = 3.5
     p_ext, _, _, _ = sample(mapping(cuts(dr_gluon=r)))
     m_gg = invariant_mass(p_ext, *GLUONS)
-    assert np.mean(m_gg < pair_floor(PT_GLUON, PT_GLUON, r)) > 0.001
-    assert np.mean(m_gg < GLUON_FLOOR) > 0.001
+    assert (m_gg < pair_floor(PT_GLUON, PT_GLUON, r)).double().mean() > 0.001
+    assert (m_gg < GLUON_FLOOR).double().mean() > 0.001
diff --git a/madspace/tests/test_rank4_kernels.py b/madspace/tests/test_rank4_kernels.py
new file mode 100644
index 000000000..e31c6a2ed
--- /dev/null
+++ b/madspace/tests/test_rank4_kernels.py
@@ -0,0 +1,22 @@
+"""An elementwise kernel on a rank-4 tensor, batch x 2 x 3 x 4.
+
+The GPU kernel launcher used to compute only the first slice of the last dimension
+of such a tensor, leaving the other elements of the output unset.
+"""
+
+import torch
+from madspace_testing import default_context, generator, randn
+
+import madspace as ms
+
+
+def test_rank4_add(n=4096):
+    rank4 = ms.Type(ms.DataType.float, ms.batch_size, [2, 3, 4])
+    fb = ms.FunctionBuilder(
+        ms.NamedTypes([("x", rank4)]), ms.NamedTypes([("y", rank4)])
+    )
+    x = fb.input(0)
+    fb.output(0, fb.add(x, x))
+    a = randn((n, 2, 3, 4), generator=generator(23))
+    (y,) = ms.FunctionRuntime(fb.function(), default_context()).call([a])
+    assert torch.equal(torch.from_dlpack(y), 2 * a)
diff --git a/madspace/tests/test_scale.py b/madspace/tests/test_scale.py
index 2548d6500..fcc6f73fb 100644
--- a/madspace/tests/test_scale.py
+++ b/madspace/tests/test_scale.py
@@ -2,9 +2,11 @@
 the dynamical scale only (as the LO `scalefact`); fixed scales are left alone.
 """
 
-import numpy as np
+import math
+
 import pytest
-from pytest import approx
+import torch
+from madspace_testing import generator, rand
 
 import madspace as ms
 
@@ -13,35 +15,35 @@
 
 def _momenta():
     """A few p p -> t t~ events, (batch, n, 4) with (E, px, py, pz)."""
-    rng = np.random.default_rng(20260917)
+    rng = generator(20260917)
     batch = 32
-    pt = rng.uniform(10.0, 400.0, size=(batch, 2))
-    phi = rng.uniform(0.0, 2 * np.pi, size=(batch, 2))
-    pz = rng.uniform(-600.0, 600.0, size=(batch, 2))
-    e = np.sqrt(M_TOP**2 + pt**2 + pz**2)
-    out = np.stack(
-        [e, pt * np.cos(phi), pt * np.sin(phi), pz], axis=-1
+    pt = 10.0 + 390.0 * rand((batch, 2), generator=rng)
+    phi = 2 * math.pi * rand((batch, 2), generator=rng)
+    pz = -600.0 + 1200.0 * rand((batch, 2), generator=rng)
+    e = torch.sqrt(M_TOP**2 + pt**2 + pz**2)
+    out = torch.stack(
+        [e, pt * torch.cos(phi), pt * torch.sin(phi), pz], dim=-1
     )  # (batch, 2, 4)
-    e_in = out[:, :, 0].sum(axis=1) / 2
-    incoming = np.zeros((batch, 2, 4))
+    e_in = out[:, :, 0].sum(dim=1) / 2
+    incoming = torch.zeros((batch, 2, 4))
     incoming[:, 0, 0] = e_in
     incoming[:, 0, 3] = e_in
     incoming[:, 1, 0] = e_in
     incoming[:, 1, 3] = -e_in
-    return np.concatenate([incoming, out], axis=1)
+    return torch.cat([incoming, out], dim=1)
 
 
 def _ht(momenta):
     """sum_i sqrt(m_i^2 + pT_i^2) over the final state."""
     final = momenta[:, 2:, :]
-    return np.sqrt(
-        np.maximum(0.0, final[..., 0] ** 2 - final[..., 3] ** 2)
-    ).sum(axis=1)
+    return torch.sqrt(
+        torch.clamp(final[..., 0] ** 2 - final[..., 3] ** 2, min=0.0)
+    ).sum(dim=1)
 
 
-def _call(scale, momenta):
-    ren, fact1, fact2 = scale(momenta)
-    return np.asarray(ren), np.asarray(fact1), np.asarray(fact2)
+def assert_approx(actual, expected, rel=1e-6, abs=1e-12):
+    """pytest.approx-like elementwise comparison of tensors."""
+    torch.testing.assert_close(actual, expected, rtol=rel, atol=abs)
 
 
 @pytest.mark.parametrize(
@@ -54,11 +56,11 @@ def _call(scale, momenta):
 def test_dynamical_scale_definitions(dyn_type, expected):
     momenta = _momenta()
     scale = ms.EnergyScale(momenta.shape[1], dyn_type, False, False, 0.0, 0.0, 0.0)
-    ren, fact1, fact2 = _call(scale, momenta)
+    ren, fact1, fact2 = scale(momenta)
     reference = expected * _ht(momenta)
-    assert ren == approx(reference)
-    assert fact1 == approx(reference)
-    assert fact2 == approx(reference)
+    assert_approx(ren, reference)
+    assert_approx(fact1, reference)
+    assert_approx(fact2, reference)
 
 
 @pytest.mark.parametrize("factor", [0.5, 1.0, 2.0])
@@ -74,11 +76,11 @@ def test_scale_factor_multiplies_the_dynamical_scale(factor):
         0.0,
         factor,
     )
-    ren, fact1, fact2 = _call(scale, momenta)
+    ren, fact1, fact2 = scale(momenta)
     reference = factor * 0.5 * _ht(momenta)
-    assert ren == approx(reference)
-    assert fact1 == approx(reference)
-    assert fact2 == approx(reference)
+    assert_approx(ren, reference)
+    assert_approx(fact1, reference)
+    assert_approx(fact2, reference)
 
 
 def test_scale_factor_default_is_one():
@@ -92,10 +94,10 @@ def test_scale_factor_default_is_one():
         0.0,
         0.0,
     )
-    without = _call(ms.EnergyScale(*args), momenta)
-    with_one = _call(ms.EnergyScale(*args, 1.0), momenta)
+    without = ms.EnergyScale(*args)(momenta)
+    with_one = ms.EnergyScale(*args, 1.0)(momenta)
     for a, b in zip(without, with_one):
-        assert a == approx(b)
+        assert_approx(a, b)
 
 
 def test_scale_factor_leaves_a_fixed_scale_alone():
@@ -111,10 +113,10 @@ def test_scale_factor_leaves_a_fixed_scale_alone():
         fixed,
         0.5,
     )
-    ren, fact1, fact2 = _call(scale, momenta)
-    assert ren == approx(np.full(momenta.shape[0], fixed))
-    assert fact1 == approx(np.full(momenta.shape[0], fixed))
-    assert fact2 == approx(np.full(momenta.shape[0], fixed))
+    ren, fact1, fact2 = scale(momenta)
+    assert_approx(ren, torch.full((momenta.shape[0],), fixed))
+    assert_approx(fact1, torch.full((momenta.shape[0],), fixed))
+    assert_approx(fact2, torch.full((momenta.shape[0],), fixed))
 
 
 def test_scale_factor_with_a_fixed_renormalisation_scale_only():
@@ -131,7 +133,7 @@ def test_scale_factor_with_a_fixed_renormalisation_scale_only():
         0.0,
         0.25,
     )
-    ren, fact1, fact2 = _call(scale, momenta)
-    assert ren == approx(np.full(momenta.shape[0], fixed))
-    assert fact1 == approx(0.25 * 0.5 * _ht(momenta))
-    assert fact2 == approx(0.25 * 0.5 * _ht(momenta))
+    ren, fact1, fact2 = scale(momenta)
+    assert_approx(ren, torch.full((momenta.shape[0],), fixed))
+    assert_approx(fact1, 0.25 * 0.5 * _ht(momenta))
+    assert_approx(fact2, 0.25 * 0.5 * _ht(momenta))
diff --git a/madspace/tests/test_sqrt_s_cut.py b/madspace/tests/test_sqrt_s_cut.py
index 4717ff6bf..805b48d33 100644
--- a/madspace/tests/test_sqrt_s_cut.py
+++ b/madspace/tests/test_sqrt_s_cut.py
@@ -26,11 +26,10 @@
 import functools
 import math
 
-import numpy as np
-import pytest
+import phasespace_helpers
+import torch
 
 import madspace as ms
-import phasespace_helpers
 
 CM_ENERGY = 13000.0
 M_TOP = 173.0
@@ -80,15 +79,15 @@ def volume(det, p):
 def test_obs_sqrt_s_is_the_partonic_invariant_mass():
     p_ext, _ = sample(mapping())
     O = ms.Observable
-    computed = np.asarray(O(PIDS, O.obs_sqrt_s, [])(p_ext)).reshape(-1)
-    assert computed == pytest.approx(s_hat(p_ext), rel=1e-9)
+    computed = O(PIDS, O.obs_sqrt_s, [])(p_ext).reshape(-1)
+    torch.testing.assert_close(computed, s_hat(p_ext), rtol=1e-9, atol=1e-12)
     # and it is the invariant mass of the final state too, i.e. it really is
     # the quantity the whole event has to clear
     out = p_ext[:, 2, :] + p_ext[:, 3, :]
-    m_out = np.sqrt(
-        np.maximum(out[:, 0] ** 2 - np.sum(out[:, 1:] ** 2, axis=1), 0.0)
+    m_out = torch.sqrt(
+        torch.clamp(out[:, 0] ** 2 - torch.sum(out[:, 1:] ** 2, dim=1), min=0.0)
     )
-    assert computed == pytest.approx(m_out, rel=1e-6)
+    torch.testing.assert_close(computed, m_out, rtol=1e-6, atol=1e-12)
 
 
 def test_the_process_has_a_threshold_of_its_own():
@@ -112,16 +111,16 @@ def test_a_cut_below_the_threshold_changes_nothing():
     the integral unchanged in the sharpest possible sense."""
     p_free, det_free = sample(mapping())
     p_cut, det_cut = sample(mapping(cuts=sqrt_s_cuts(NO_OP_CUT)))
-    assert np.array_equal(det_free, det_cut)
-    assert np.array_equal(p_free, p_cut)
+    assert torch.equal(det_free, det_cut)
+    assert torch.equal(p_free, p_cut)
 
 
 def test_a_cut_at_zero_changes_nothing():
     """The same statement for the value the run card ships with."""
     p_free, det_free = sample(mapping())
     p_cut, det_cut = sample(mapping(cuts=sqrt_s_cuts(0.0)))
-    assert np.array_equal(det_free, det_cut)
-    assert np.array_equal(p_free, p_cut)
+    assert torch.equal(det_free, det_cut)
+    assert torch.equal(p_free, p_cut)
 
 
 # --------------------------------------------------------------------------
@@ -138,7 +137,7 @@ def test_cut_bounds_the_sampled_s_hat():
     # the region has not collapsed onto its own boundary
     assert s.max() > 2 * BITING_CUT
     # and the cut is no longer filtering anything out
-    assert np.all(det != 0.0)
+    assert torch.all(det != 0.0)
 
 
 def test_without_the_bound_the_same_region_is_wasted():
@@ -147,7 +146,7 @@ def test_without_the_bound_the_same_region_is_wasted():
     what the boundary removes, so if it ever stops being true the test above
     has stopped proving anything."""
     p_ext, _ = sample(mapping())
-    assert np.mean(s_hat(p_ext) < BITING_CUT) > 0.25
+    assert (s_hat(p_ext) < BITING_CUT).double().mean() > 0.25
 
 
 def test_cut_volume_is_unchanged_by_the_boundary():
@@ -156,13 +155,14 @@ def test_cut_volume_is_unchanged_by_the_boundary():
     cut or only the filtering does."""
     p_free, det_free = sample(mapping(), seed=SEED)
     passes = s_hat(p_free) >= BITING_CUT
-    external = np.where(passes, volume(det_free, p_free), 0.0)
+    external = torch.where(passes, volume(det_free, p_free), 0.0)
 
     p_cut, det_cut = sample(mapping(cuts=sqrt_s_cuts(BITING_CUT)), seed=SEED + 1)
     piped = volume(det_cut, p_cut)
 
     error = math.sqrt(
-        external.std() ** 2 / len(external) + piped.std() ** 2 / len(piped)
+        external.std(correction=0).item() ** 2 / len(external)
+        + piped.std(correction=0).item() ** 2 / len(piped)
     )
     assert error > 0.0
     assert abs(external.mean() - piped.mean()) < 5.0 * error
@@ -175,8 +175,8 @@ def test_the_cut_removes_more_than_the_comparison_could_hide():
     there, not merely nonzero."""
     p_free, det_free = sample(mapping())
     total = volume(det_free, p_free)
-    removed = np.where(s_hat(p_free) < BITING_CUT, total, 0.0)
-    error = total.std() / math.sqrt(len(total))
+    removed = torch.where(s_hat(p_free) < BITING_CUT, total, 0.0)
+    error = total.std(correction=0).item() / math.sqrt(len(total))
     assert error > 0.0
     assert removed.mean() > 10.0 * error
 
@@ -190,9 +190,9 @@ def test_a_cut_above_the_beam_energy_is_empty_and_finite():
     """Such a cut leaves nothing to sample at all. It must come out as an empty
     result, not as an inverted sampling range with nan weights."""
     p_ext, det = sample(mapping(cuts=sqrt_s_cuts(2 * CM_ENERGY)))
-    assert np.all(det == 0.0)
-    assert np.all(np.isfinite(det))
-    assert np.all(np.isfinite(p_ext))
+    assert torch.all(det == 0.0)
+    assert torch.all(torch.isfinite(det))
+    assert torch.all(torch.isfinite(p_ext))
 
 
 def test_leptonic_s_hat_is_fixed_so_the_cut_only_filters():
@@ -200,8 +200,9 @@ def test_leptonic_s_hat_is_fixed_so_the_cut_only_filters():
     luminosity to sample and hence no range to narrow, so the cut stays a plain
     filter and must still be applied as one."""
     p_ext, det = sample(mapping(leptonic=True))
-    assert s_hat(p_ext) == pytest.approx(CM_ENERGY, rel=1e-9)
+    s = s_hat(p_ext)
+    torch.testing.assert_close(s, torch.full_like(s, CM_ENERGY), rtol=1e-9, atol=1e-12)
     _, det_kept = sample(mapping(cuts=sqrt_s_cuts(BITING_CUT), leptonic=True))
-    assert np.array_equal(det_kept, det)
+    assert torch.equal(det_kept, det)
     _, det_gone = sample(mapping(cuts=sqrt_s_cuts(2 * CM_ENERGY), leptonic=True))
-    assert np.all(det_gone == 0.0)
+    assert torch.all(det_gone == 0.0)
diff --git a/madspace/tests/test_survey_job_failure.py b/madspace/tests/test_survey_job_failure.py
index d56a6816f..998d9dd95 100644
--- a/madspace/tests/test_survey_job_failure.py
+++ b/madspace/tests/test_survey_job_failure.py
@@ -18,8 +18,13 @@
 import sys
 import textwrap
 
+import pytest
+
 import madspace as ms
 
+# the survey runs on CPU threads of its own context
+pytestmark = pytest.mark.cpu_only
+
 TIMEOUT = 120
 
 SCRIPT = textwrap.dedent(
diff --git a/madspace/tests/test_systematics.py b/madspace/tests/test_systematics.py
index 31732a0c6..bdc49e247 100644
--- a/madspace/tests/test_systematics.py
+++ b/madspace/tests/test_systematics.py
@@ -38,9 +38,13 @@ def _pdf_dir():
 
 
 PDF_DIR = _pdf_dir()
-pytestmark = pytest.mark.skipif(
-    PDF_DIR is None, reason=f"pdf set {PDF_SET} required to run this test"
-)
+# host-side: LHAPDF, the systematics calculator and explicit CPU contexts
+pytestmark = [
+    pytest.mark.skipif(
+        PDF_DIR is None, reason=f"pdf set {PDF_SET} required to run this test"
+    ),
+    pytest.mark.cpu_only,
+]
 
 
 def grid_files(set_name=PDF_SET, member=0):
@@ -88,9 +92,7 @@ def ref(grids):
 def sample_points(grid, n, seed=1):
     rng = np.random.default_rng(seed)
     x = 10 ** rng.uniform(np.log10(grid.x[0]) + 1e-6, 0, n)
-    q = 10 ** rng.uniform(
-        np.log10(grid.q[0]) + 1e-6, np.log10(grid.q[-1]) - 1e-6, n
-    )
+    q = 10 ** rng.uniform(np.log10(grid.q[0]) + 1e-6, np.log10(grid.q[-1]) - 1e-6, n)
     return x, q
 
 
@@ -159,8 +161,14 @@ def test_scale_variations(grids, ref):
     variations = calc.variations
     assert len(variations) == 8
     assert [(v.mur, v.muf) for v in variations] == [
-        (0.5, 0.5), (0.5, 1.0), (0.5, 2.0), (1.0, 0.5),
-        (1.0, 2.0), (2.0, 0.5), (2.0, 1.0), (2.0, 2.0),
+        (0.5, 0.5),
+        (0.5, 1.0),
+        (0.5, 2.0),
+        (1.0, 0.5),
+        (1.0, 2.0),
+        (2.0, 0.5),
+        (2.0, 1.0),
+        (2.0, 2.0),
     ]
     assert calc.weight_ids == list(range(1, 9))
     events = make_events(ref)
@@ -191,7 +199,10 @@ def test_identities(grids, ref):
     # leptonic beams: mu_F variations are dropped, mu_R ones remain
     config = make_config(has_pdf=False, pdf_members=[])
     calc = ms.SystematicsCalculator(
-        config, [ms.SubprocessSystArgs(qcd_power=2, beam_pdgs=[[11, -11]])], None, alpha_s
+        config,
+        [ms.SubprocessSystArgs(qcd_power=2, beam_pdgs=[[11, -11]])],
+        None,
+        alpha_s,
     )
     assert [(v.mur, v.muf) for v in calc.variations] == [(0.5, 1.0), (2.0, 1.0)]
 
@@ -200,7 +211,10 @@ def test_non_uniform_qcd_power_drops_mur(grids):
     grid, alpha_s = grids
     config = make_config()
     calc = ms.SystematicsCalculator(
-        config, [ms.SubprocessSystArgs(qcd_power=-1, beam_pdgs=[[2, -2]])], grid, alpha_s
+        config,
+        [ms.SubprocessSystArgs(qcd_power=-1, beam_pdgs=[[2, -2]])],
+        grid,
+        alpha_s,
     )
     assert [(v.mur, v.muf) for v in calc.variations] == [(1.0, 0.5), (1.0, 2.0)]
     assert any("renormalisation scale" in w for w in calc.warnings)
@@ -213,12 +227,22 @@ def test_pdf_member_and_header(grids, ref):
     grid, alpha_s = grids
     grid_file, info_file = grid_files()
     other = ms.PdfMemberSpec(
-        set_name="OtherSet", set_lhaid=900000, member=1, grid_file=grid_file,
-        info_file=info_file, error_type="hessian", description="a => b > c",
+        set_name="OtherSet",
+        set_lhaid=900000,
+        member=1,
+        grid_file=grid_file,
+        info_file=info_file,
+        error_type="hessian",
+        description="a => b > c",
     )
     same = ms.PdfMemberSpec(
-        set_name=PDF_SET, set_lhaid=338500, member=1, grid_file=grid_file,
-        info_file=info_file, error_type="replicas", description="",
+        set_name=PDF_SET,
+        set_lhaid=338500,
+        member=1,
+        grid_file=grid_file,
+        info_file=info_file,
+        error_type="replicas",
+        description="",
     )
     config = make_config(mur=[1.0], muf=[1.0], pdf_members=[same, other])
     calc = ms.SystematicsCalculator(
@@ -232,10 +256,19 @@ def test_pdf_member_and_header(grids, ref):
         assert row == approx([w0, w0, w0], rel=1e-12)
     header = calc.initrwgt()
     assert '<weightgroup name="%s" combine="replicas">' % PDF_SET in header
-    assert '<weightgroup name="OtherSet" combine="hessian"> # 900000: a ; b .gt. c' in header
+    assert (
+        '<weightgroup name="OtherSet" combine="hessian"> # 900000: a ; b .gt. c'
+        in header
+    )
     assert '<weight id="1" MUR="1.0" MUF="1.0" PDF="338500" >  </weight>' in header
-    assert '<weight id="2" MUR="1.0" MUF="1.0" PDF="338501" > PDF=338500 MemberID=1 </weight>' in header
-    assert '<weight id="3" MUR="1.0" MUF="1.0" PDF="900001" > PDF=900000 MemberID=1 </weight>' in header
+    assert (
+        '<weight id="2" MUR="1.0" MUF="1.0" PDF="338501" > PDF=338500 MemberID=1 </weight>'
+        in header
+    )
+    assert (
+        '<weight id="3" MUR="1.0" MUF="1.0" PDF="900001" > PDF=900000 MemberID=1 </weight>'
+        in header
+    )
     assert header.count("<weightgroup") == header.count("</weightgroup>") == 2
     summary = json.loads(calc.summary())
     assert [v["id"] for v in summary["variations"]] == [1, 2, 3]
@@ -248,7 +281,10 @@ def test_lhe_event_rwgt_format():
     event.rwgt_ids = [1, 2]
     event.rwgt = [1.5e3, -2.0]
     text = event.format()
-    assert "<rwgt>\n<wgt id='1'> +1.5000000e+03 </wgt>\n<wgt id='2'> -2.0000000e+00 </wgt>\n</rwgt>\n</event>" in text
+    assert (
+        "<rwgt>\n<wgt id='1'> +1.5000000e+03 </wgt>\n<wgt id='2'> -2.0000000e+00 </wgt>\n</rwgt>\n</event>"
+        in text
+    )
 
 
 def test_config_json_roundtrip():
@@ -278,11 +314,17 @@ def random_momenta(n, seed=5):
         s = 4 * e1 * e2
         pcm = np.sqrt(s) / 2
         theta, phi = np.arccos(rng.uniform(-1, 1)), rng.uniform(0, 2 * np.pi)
-        px, py, pz = pcm * np.sin(theta) * np.cos(phi), pcm * np.sin(theta) * np.sin(phi), pcm * np.cos(theta)
+        px, py, pz = (
+            pcm * np.sin(theta) * np.cos(phi),
+            pcm * np.sin(theta) * np.sin(phi),
+            pcm * np.cos(theta),
+        )
         beta = (e1 - e2) / (e1 + e2)
         gamma = 1 / np.sqrt(1 - beta**2)
+
         def boost(e, pz):
             return gamma * (e + beta * pz), gamma * (pz + beta * e)
+
         ea, pza = boost(pcm, pz)
         eb, pzb = boost(pcm, -pz)
         out.append(p_in + [[ea, px, py, pza], [eb, -px, -py, pzb]])
@@ -297,13 +339,24 @@ def test_dynamical_scales(grids, ref):
     for p in momenta[:5]:
         p = np.array(p)
         fin = p[2:]
-        et = sum(e * np.sqrt((px**2 + py**2) / (px**2 + py**2 + pz**2)) for e, px, py, pz in fin)
+        et = sum(
+            e * np.sqrt((px**2 + py**2) / (px**2 + py**2 + pz**2))
+            for e, px, py, pz in fin
+        )
         ht = sum(np.sqrt(max(0.0, e**2 - pz**2)) for e, px, py, pz in fin)
         shat = np.sqrt((p[0, 0] + p[1, 0]) ** 2 - (p[0, 3] + p[1, 3]) ** 2)
-        assert ms.SystematicsCalculator.dynamical_scale(1, p.tolist()) == approx(et, rel=1e-12)
-        assert ms.SystematicsCalculator.dynamical_scale(2, p.tolist()) == approx(ht, rel=1e-12)
-        assert ms.SystematicsCalculator.dynamical_scale(3, p.tolist()) == approx(ht / 2, rel=1e-12)
-        assert ms.SystematicsCalculator.dynamical_scale(4, p.tolist()) == approx(shat, rel=1e-12)
+        assert ms.SystematicsCalculator.dynamical_scale(1, p.tolist()) == approx(
+            et, rel=1e-12
+        )
+        assert ms.SystematicsCalculator.dynamical_scale(2, p.tolist()) == approx(
+            ht, rel=1e-12
+        )
+        assert ms.SystematicsCalculator.dynamical_scale(3, p.tolist()) == approx(
+            ht / 2, rel=1e-12
+        )
+        assert ms.SystematicsCalculator.dynamical_scale(4, p.tolist()) == approx(
+            shat, rel=1e-12
+        )
 
     config = make_config(mur=[0.5, 1.0], muf=[1.0, 2.0], dyn_scales=[3, 4])
     calc = ms.SystematicsCalculator(
@@ -312,9 +365,17 @@ def test_dynamical_scales(grids, ref):
     # (mur, muf) grid at the generated scale minus the nominal point, then the
     # full grid for each dynamical choice
     assert [(v.mur, v.muf, v.dyn) for v in calc.variations] == [
-        (0.5, 1.0, -1), (0.5, 2.0, -1), (1.0, 2.0, -1),
-        (0.5, 1.0, 3), (0.5, 2.0, 3), (1.0, 1.0, 3), (1.0, 2.0, 3),
-        (0.5, 1.0, 4), (0.5, 2.0, 4), (1.0, 1.0, 4), (1.0, 2.0, 4),
+        (0.5, 1.0, -1),
+        (0.5, 2.0, -1),
+        (1.0, 2.0, -1),
+        (0.5, 1.0, 3),
+        (0.5, 2.0, 3),
+        (1.0, 1.0, 3),
+        (1.0, 2.0, 3),
+        (0.5, 1.0, 4),
+        (0.5, 2.0, 4),
+        (1.0, 1.0, 4),
+        (1.0, 2.0, 4),
     ]
     header = calc.initrwgt()
     assert 'DYN_SCALE="3" PDF="338500" > dyn_scale_choice=HT/2' in header
@@ -326,9 +387,17 @@ def test_dynamical_scales(grids, ref):
         x1, x2 = events["x1"][i], events["x2"][i]
         nominal_pdf = ref.pdf(2, x1, mu) * ref.pdf(2, x2, mu)
         for var, w in zip(calc.variations, weights[i]):
-            mu_dyn = mu if var.dyn == -1 else ms.SystematicsCalculator.dynamical_scale(var.dyn, momenta[i])
+            mu_dyn = (
+                mu
+                if var.dyn == -1
+                else ms.SystematicsCalculator.dynamical_scale(var.dyn, momenta[i])
+            )
             r_alpha = (ref.alphas(var.mur * mu_dyn) / ref.alphas(mu)) ** 2
-            r_pdf = ref.pdf(2, x1, var.muf * mu_dyn) * ref.pdf(2, x2, var.muf * mu_dyn) / nominal_pdf
+            r_pdf = (
+                ref.pdf(2, x1, var.muf * mu_dyn)
+                * ref.pdf(2, x2, var.muf * mu_dyn)
+                / nominal_pdf
+            )
             assert w == approx(w0 * r_alpha * r_pdf, rel=1e-12)
 
 
@@ -341,17 +410,29 @@ def test_event_histograms(grids, ref):
     pt = np.array([np.hypot(p[2][1], p[2][2]) for p in momenta])
     ctx = ms.Context(device=ms.cpu_device(), thread_count=1)
     pids = [2, 2, 2, 2]
-    obs = ms.ObservableValues([
-        ms.Observable(pids, observable="pt", select_pids=[[2]], order_observable="pt",
-                      order_indices=[1], name="jet-pt"),
-        ms.Observable(pids, observable="sqrt_s", select_pids=[], name="sqrt_s"),
-    ])
-    specs = [ms.EventHistogramSpec("jet-pt", 0.0, 500.0, 10),
-             ms.EventHistogramSpec("sqrt_s", 0.0, 2000.0, 8)]
+    obs = ms.ObservableValues(
+        [
+            ms.Observable(
+                pids,
+                observable="pt",
+                select_pids=[[2]],
+                order_observable="pt",
+                order_indices=[1],
+                name="jet-pt",
+            ),
+            ms.Observable(pids, observable="sqrt_s", select_pids=[], name="sqrt_s"),
+        ]
+    )
+    specs = [
+        ms.EventHistogramSpec("jet-pt", 0.0, 500.0, 10),
+        ms.EventHistogramSpec("sqrt_s", 0.0, 2000.0, 8),
+    ]
     hists = ms.EventHistograms(ctx, specs, [ms.SubprocessObservables(obs, 4)])
     rng = np.random.default_rng(7)
     weight = np.full(n, 3.0)
-    syst = np.stack([weight * 0.8, weight * 1.3, weight], axis=1)  # two scale-like columns + nominal copy
+    syst = np.stack(
+        [weight * 0.8, weight * 1.3, weight], axis=1
+    )  # two scale-like columns + nominal copy
     hists.fill(list(weight), [0] * n, momenta, syst.tolist())
     data = json.loads(hists.to_json())
     assert [h["name"] for h in data] == ["jet-pt", "sqrt_s"]
@@ -363,7 +444,9 @@ def test_event_histograms(grids, ref):
     assert jet["bin_values"][1:-1] == approx(list(3.0 * counts / n))
     # binomial error of an unweighted sample
     m = counts[counts.argmax()]
-    assert jet["bin_errors"][1 + counts.argmax()] == approx(3.0 * np.sqrt(m * (1 - m / n)) / n)
+    assert jet["bin_errors"][1 + counts.argmax()] == approx(
+        3.0 * np.sqrt(m * (1 - m / n)) / n
+    )
     # a systematics calculator adds the envelope: fake one with two scale variations
     config = make_config(mur=[0.5, 1.0, 2.0], muf=[1.0], pdf_members=[])
     calc = ms.SystematicsCalculator(
diff --git a/madspace/tests/test_t_channel.py b/madspace/tests/test_t_channel.py
index 1040fd560..6e6d18e32 100644
--- a/madspace/tests/test_t_channel.py
+++ b/madspace/tests/test_t_channel.py
@@ -1,12 +1,18 @@
 import math
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand, randint
 from pytest import approx
 
 import madspace as ms
 
-np.set_printoptions(linewidth=1000)
+torch.set_printoptions(linewidth=1000)
+
+
+def assert_approx(actual, expected, rel=1e-6, abs=1e-12):
+    """pytest.approx-like elementwise comparison of tensors."""
+    torch.testing.assert_close(actual, expected, rtol=rel, atol=abs)
 
 
 def _color_order(n_out, mode):
@@ -18,19 +24,19 @@ def _color_order(n_out, mode):
 
 def _fwd_inputs(mapping, rng, n):
     """[random] (+ [discrete] when the mode declares discrete choices)."""
-    r = rng.random((n, mapping.random_dim()))
+    r = rand((n, mapping.random_dim()), generator=rng)
     inputs = [r]
     disc = None
     dd = mapping.discrete_dim()
     if dd:
-        disc = rng.integers(0, 2, size=(n, dd)).astype(np.int32)
+        disc = randint(0, 2, (n, dd), generator=rng, dtype=torch.int32)
         inputs.append(disc)
     return inputs, r, disc
 
 
 @pytest.fixture
 def rng():
-    return np.random.default_rng(1234)
+    return generator(43)
 
 
 @pytest.fixture(
@@ -100,9 +106,9 @@ def test_t_channel_masses(masses, rng, mode):
         p_ext = p_ext[physical_mask]
         batch_phys = p_ext.shape[0]
 
-    m_ext_true = np.full((batch_phys, len(masses)), masses)
-    m_ext = np.sqrt(
-        np.maximum(0, p_ext[:, :, 0] ** 2 - np.sum(p_ext[:, :, 1:] ** 2, axis=2))
+    m_ext_true = torch.tensor(masses).expand(batch_phys, -1)
+    m_ext = torch.sqrt(
+        torch.clamp(p_ext[:, :, 0] ** 2 - torch.sum(p_ext[:, :, 1:] ** 2, dim=2), min=0)
     )
     # The massless on-shell reconstruction m = sqrt(E^2 - |p|^2) has an FP floor
     # that scales with the particle energy (~2e-6 * E); the deep color_ordered
@@ -111,8 +117,10 @@ def test_t_channel_masses(masses, rng, mode):
     # energy-scaled term (a genuine kinematic error would exceed it ~1000x).
     abs_tol = 1e-3
     if mode == ms.PhaseSpaceMapping.color_ordered:
-        abs_tol = 1e-3 + 2e-6 * np.abs(p_ext[:, :, 0])
-    assert np.all(np.abs(m_ext - m_ext_true) <= abs_tol + 1e-3 * np.abs(m_ext_true))
+        abs_tol = 1e-3 + 2e-6 * torch.abs(p_ext[:, :, 0])
+    assert torch.all(
+        torch.abs(m_ext - m_ext_true) <= abs_tol + 1e-3 * torch.abs(m_ext_true)
+    )
 
 
 def test_t_channel_incoming(masses, rng, mode):
@@ -133,18 +141,21 @@ def test_t_channel_incoming(masses, rng, mode):
         p_ext = p_ext[physical_mask]
         batch_phys = p_ext.shape[0]
 
-    zeros = np.zeros(batch_phys)
+    zeros = torch.zeros(batch_phys)
     p_a = p_ext[:, 0]
     p_b = p_ext[:, 1]
     e_beam = 0.5 * CM_ENERGY
 
-    assert p_a[:, 0] == approx(p_a[:, 3]) and p_b[:, 0] == approx(-p_b[:, 3])
-    assert p_a[:, 1] == approx(zeros) and p_a[:, 2] == approx(zeros)
-    assert p_b[:, 1] == approx(zeros) and p_b[:, 2] == approx(zeros)
-    assert np.all(x1 >= 0) and np.all(x1 <= 1)
-    assert np.all(x2 >= 0) and np.all(x2 <= 1)
-    assert p_a[:, 0] == approx(e_beam * x1)
-    assert p_b[:, 0] == approx(e_beam * x2)
+    assert_approx(p_a[:, 0], p_a[:, 3])
+    assert_approx(p_b[:, 0], -p_b[:, 3])
+    assert_approx(p_a[:, 1], zeros)
+    assert_approx(p_a[:, 2], zeros)
+    assert_approx(p_b[:, 1], zeros)
+    assert_approx(p_b[:, 2], zeros)
+    assert torch.all(x1 >= 0) and torch.all(x1 <= 1)
+    assert torch.all(x2 >= 0) and torch.all(x2 <= 1)
+    assert_approx(p_a[:, 0], e_beam * x1)
+    assert_approx(p_b[:, 0], e_beam * x2)
 
 
 def test_t_channel_momentum_conservation(masses, rng, mode):
@@ -161,10 +172,11 @@ def test_t_channel_momentum_conservation(masses, rng, mode):
         physical_mask = det != 0.0
         p_ext = p_ext[physical_mask]
 
-    p_in = np.sum(p_ext[:, :2], axis=1)
-    p_out = np.sum(p_ext[:, 2:], axis=1)
+    p_in = torch.sum(p_ext[:, :2], dim=1)
+    p_out = torch.sum(p_ext[:, 2:], dim=1)
 
-    assert p_out == approx(p_in, rel=1e-6, abs=1e-9)
+    # the transverse components are 0 up to the rounding of momenta of order CM_ENERGY
+    assert_approx(p_out, p_in, rel=1e-6, abs=1e-12 * CM_ENERGY)
 
 
 def test_t_channel_inverse(masses, rng, mode):
@@ -191,14 +203,13 @@ def test_t_channel_inverse(masses, rng, mode):
     out = mapping.map_inverse((p_ext, x1, x2))
     det_inv = out[-1]
     r_inv = out[0]
-    one_batch = np.ones_like(det)
+    one_batch = torch.ones_like(det)
     # Continuous randoms round-trip; the discrete int channel (color_ordered)
     # is checked for exact bin/solution recovery separately.
-    assert r_inv == approx(r, abs=1e-3, rel=1e-3)
-    assert det * det_inv == approx(one_batch, rel=1e-5)
+    assert_approx(r_inv, r, abs=1e-3, rel=1e-3)
+    assert_approx(det * det_inv, one_batch, rel=1e-5)
     if disc is not None:
-        disc_inv = np.asarray(out[1]).astype(np.int64)
-        assert np.array_equal(disc_inv, disc.astype(np.int64))
+        assert torch.equal(out[1].long(), disc.long())
 
 
 @pytest.mark.parametrize(
@@ -241,8 +252,8 @@ def test_t_channel_phase_space_volume(particle_count, energy, rng, mode):
     )
     if mode == ms.PhaseSpaceMapping.chili:
         ps_volume /= (particle_count - 1) ** 2  # integration over x1*x2 in [0,1]
-    std_error = np.std(det) / np.sqrt(sample_count)
-    assert np.mean(det) == approx(ps_volume, abs=3 * std_error, rel=1e-6)
+    std_error = det.std(correction=0).item() / math.sqrt(sample_count)
+    assert det.mean().item() == approx(ps_volume, abs=3 * std_error, rel=1e-6)
 
 
 # ---------------------------------------------------------------------------
@@ -288,26 +299,28 @@ def test_color_ordered_color_orders(rng, color_order, out_masses, _label):
     p_ext, x1, x2, det = mapping.map_forward(inputs)
 
     # (1) momentum conservation
-    p_in = np.sum(p_ext[:, :2], axis=1)
-    p_out = np.sum(p_ext[:, 2:], axis=1)
-    assert p_out == approx(p_in, rel=1e-6, abs=1e-9)
+    p_in = torch.sum(p_ext[:, :2], dim=1)
+    p_out = torch.sum(p_ext[:, 2:], dim=1)
+    # the transverse components are 0 up to the rounding of momenta of order CM_ENERGY
+    assert_approx(p_out, p_in, rel=1e-6, abs=1e-12 * CM_ENERGY)
 
     # (2) on-shell external masses (energy-scaled FP floor for deep boost chains)
-    m_ext_true = np.full((BATCH_SIZE, len(masses)), masses)
-    m_ext = np.sqrt(
-        np.maximum(0, p_ext[:, :, 0] ** 2 - np.sum(p_ext[:, :, 1:] ** 2, axis=2))
+    m_ext_true = torch.tensor(masses).expand(BATCH_SIZE, -1)
+    m_ext = torch.sqrt(
+        torch.clamp(p_ext[:, :, 0] ** 2 - torch.sum(p_ext[:, :, 1:] ** 2, dim=2), min=0)
+    )
+    abs_tol = 1e-3 + 2e-6 * torch.abs(p_ext[:, :, 0])
+    assert torch.all(
+        torch.abs(m_ext - m_ext_true) <= abs_tol + 1e-3 * torch.abs(m_ext_true)
     )
-    abs_tol = 1e-3 + 2e-6 * np.abs(p_ext[:, :, 0])
-    assert np.all(np.abs(m_ext - m_ext_true) <= abs_tol + 1e-3 * np.abs(m_ext_true))
 
     # (3) inverse: continuous randoms round-trip, discrete choices exactly,
     out = mapping.map_inverse((p_ext, x1, x2))
     r_inv, det_inv = out[0], out[-1]
-    assert r_inv == approx(r, abs=2e-3, rel=1e-3)
-    rt_err = np.abs(det * det_inv - 1.0)
-    assert np.quantile(rt_err, 0.99) < 1e-5
-    assert np.isfinite(rt_err).all()
-    assert np.mean(rt_err > 1e-2) < 0.01
+    assert_approx(r_inv, r, abs=2e-3, rel=1e-3)
+    rt_err = torch.abs(det * det_inv - 1.0)
+    assert torch.quantile(rt_err, 0.99) < 1e-5
+    assert torch.isfinite(rt_err).all()
+    assert torch.mean((rt_err > 1e-2).double()) < 0.01
     if disc is not None:
-        disc_inv = np.asarray(out[1]).astype(np.int64)
-        assert np.array_equal(disc_inv, disc.astype(np.int64))
+        assert torch.equal(out[1].long(), disc.long())
diff --git a/madspace/tests/test_three_body_decay.py b/madspace/tests/test_three_body_decay.py
index dec518f8a..bede8287b 100644
--- a/madspace/tests/test_three_body_decay.py
+++ b/madspace/tests/test_three_body_decay.py
@@ -1,15 +1,21 @@
+import math
 from collections import namedtuple
 
-import numpy as np
 import pytest
-from pytest import approx
+import torch
+from madspace_testing import generator, rand, randn
 
 import madspace as ms
 
 
+def assert_approx(actual, expected):
+    """Elementwise pytest.approx with its default tolerances, for tensors."""
+    torch.testing.assert_close(actual, expected, rtol=1e-6, atol=1e-12)
+
+
 @pytest.fixture
 def rng():
-    return np.random.default_rng(1234)
+    return generator(1234)
 
 
 @pytest.fixture(
@@ -24,13 +30,13 @@ def rng():
 )
 def mapping_and_args(request):
     com = request.param["com"]
-    zeros = np.zeros(N)
+    zeros = torch.zeros(N)
     if request.param["decay"]:
         mapping = ms.ThreeBodyDecay(com=com)
         if com:
 
             def make_args(point):
-                p0 = np.stack([point.m0, zeros, zeros, zeros], axis=1)
+                p0 = torch.stack([point.m0, zeros, zeros, zeros], dim=1)
                 return (
                     [
                         point.r1,
@@ -73,7 +79,7 @@ def make_args(point):
 
 
 def mass(momentum):
-    return np.sqrt(momentum[:, 0] ** 2 - np.sum(momentum[:, 1:] ** 2, axis=1))
+    return torch.sqrt(momentum[:, 0] ** 2 - torch.sum(momentum[:, 1:] ** 2, dim=1))
 
 
 InputPoint = namedtuple(
@@ -84,31 +90,35 @@ def mass(momentum):
 N = 10000
 
 
+def uniform(rng, low, high):
+    return low + (high - low) * rand(N, generator=rng)
+
+
 @pytest.fixture
 def input_points(rng):
     max_mass = 125.0
-    r1 = rng.random(N)
-    r2 = rng.random(N)
-    r3 = rng.random(N)
-    r4 = rng.random(N)
-    r5 = rng.random(N)
-
-    pa = np.array([[0, 0, 4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-    ma = rng.uniform(2.5 * max_mass, 1000.0, N)
-    ea = np.sqrt(np.sum(pa**2, axis=1) + ma**2)
-    pa = np.concatenate([ea[:, None], pa], axis=1)
-
-    pb = np.array([[0, 0, -4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-    mb = rng.uniform(20.0, 500.0, N)
-    eb = np.sqrt(np.sum(pb**2, axis=1) + mb**2)
-    pb = np.concatenate([eb[:, None], pb], axis=1)
+    r1 = rand(N, generator=rng)
+    r2 = rand(N, generator=rng)
+    r3 = rand(N, generator=rng)
+    r4 = rand(N, generator=rng)
+    r5 = rand(N, generator=rng)
+
+    pa = torch.tensor([[0, 0, 4000.0]]) + 500.0 * randn((N, 3), generator=rng)
+    ma = uniform(rng, 2.5 * max_mass, 1000.0)
+    ea = torch.sqrt(torch.sum(pa**2, dim=1) + ma**2)
+    pa = torch.cat([ea[:, None], pa], dim=1)
+
+    pb = torch.tensor([[0, 0, -4000.0]]) + 500.0 * randn((N, 3), generator=rng)
+    mb = uniform(rng, 20.0, 500.0)
+    eb = torch.sqrt(torch.sum(pb**2, dim=1) + mb**2)
+    pb = torch.cat([eb[:, None], pb], dim=1)
 
     p0 = pa + pb
     m0 = mass(p0)
 
-    m1 = rng.uniform(1.0, max_mass, N)
-    m2 = rng.uniform(1.0, max_mass, N)
-    m3 = rng.uniform(1.0, max_mass, N)
+    m1 = uniform(rng, 1.0, max_mass)
+    m2 = uniform(rng, 1.0, max_mass)
+    m3 = uniform(rng, 1.0, max_mass)
 
     return InputPoint(r1, r2, r3, r4, r5, m0, m1, m2, m3, p0, pa, pb)
 
@@ -117,7 +127,7 @@ def test_momentum_conservation(mapping_and_args, input_points):
     mapping, make_args, com = mapping_and_args
     inputs, conditions, p0 = make_args(input_points)
     p1, p2, p3, det = mapping.map_forward(inputs, conditions)
-    assert p1 + p2 + p3 == approx(p0)
+    assert_approx(p1 + p2 + p3, p0)
 
 
 def test_inverse(mapping_and_args, input_points):
@@ -125,9 +135,11 @@ def test_inverse(mapping_and_args, input_points):
     inputs, conditions, p0 = make_args(input_points)
     p1, p2, p3, det = mapping.map_forward(inputs, conditions)
     *inv_inputs, inv_det = mapping.map_inverse([p1, p2, p3], conditions)
-    assert inv_det == approx(1 / det)
+    assert_approx(inv_det, 1 / det)
     for i, (inp, inv_inp) in enumerate(zip(inputs, inv_inputs)):
-        assert inp == approx(inv_inp), f"Mismatch at index {i}: {inp} vs {inv_inp}"
+        torch.testing.assert_close(
+            inp, inv_inp, rtol=1e-6, atol=1e-12, msg=f"Mismatch at index {i}"
+        )
 
 
 def test_outgoing_masses(mapping_and_args, input_points):
@@ -141,39 +153,39 @@ def test_outgoing_masses(mapping_and_args, input_points):
     index = 5
     print(p3[index])
     print(m3[index], input_points.m3[index])
-    assert m0 == approx(input_points.m0)
-    assert m1 == approx(input_points.m1)
-    assert m2 == approx(input_points.m2)
-    assert m3 == approx(input_points.m3)
+    assert_approx(m0, input_points.m0)
+    assert_approx(m1, input_points.m1)
+    assert_approx(m2, input_points.m2)
+    assert_approx(m3, input_points.m3)
 
 
 def test_phase_space_via_two_2body(mapping_and_args, rng):
     mapping_13, make_args, com = mapping_and_args  # your ThreeBodyDecay
     N = 50000
-    r = rng.random((N, 5))  # 5 random numbers for 1->3 decay
+    r = rand((N, 5), generator=rng)  # 5 random numbers for 1->3 decay
     u = r[:, 0]  # first random number for s12
 
     # Fix parent in COM for clarity (or use your fixture’s p0)
     m0, m1, m2, m3 = 120.0, 10.0, 15.0, 5.0
-    zeros = np.zeros(N)
-    p0 = np.stack([np.full(N, m0), zeros, zeros, zeros], axis=1)
+    zeros = torch.zeros(N)
+    p0 = torch.stack([torch.full((N,), m0), zeros, zeros, zeros], dim=1)
 
     smin = (m1 + m2) ** 2
     smax = (m0 - m3) ** 2
     s12 = smin + u * (smax - smin)
-    mQ = np.sqrt(s12)
+    mQ = torch.sqrt(s12)
 
     # --- First 2-body: 0 -> Q + 3 ---
     mapping1 = ms.TwoBodyDecay(com=True)
-    inputs1 = [r[:, 1], r[:, 2], np.full(N, m0), np.full(N, mQ), np.full(N, m3)]
+    inputs1 = [r[:, 1], r[:, 2], torch.full((N,), m0), mQ, torch.full((N,), m3)]
     pQ, p3, det1 = mapping1.map_forward(inputs1, [])
 
     # Ensure pQ has invariant mass sqrt(s12)
-    assert approx(np.sqrt(pQ[:, 0] ** 2 - (pQ[:, 1:] ** 2).sum(1))) == mQ
+    assert_approx(torch.sqrt(pQ[:, 0] ** 2 - (pQ[:, 1:] ** 2).sum(1)), mQ)
 
     # --- Second 2-body: Q -> 1 + 2 ---
     mapping2 = ms.TwoBodyDecay(com=False)  # use actual pQ (boosted)
-    inputs2 = [r[:, 3], r[:, 4], mQ, np.full(N, m1), np.full(N, m2), pQ]
+    inputs2 = [r[:, 3], r[:, 4], mQ, torch.full((N,), m1), torch.full((N,), m2), pQ]
     p1, p2, det2 = mapping2.map_forward(inputs2, [])
 
     # Reference Jacobian
@@ -187,10 +199,10 @@ def test_phase_space_via_two_2body(mapping_and_args, rng):
             r[:, 2],
             r[:, 3],
             r[:, 4],
-            np.full(N, m0),
-            np.full(N, m1),
-            np.full(N, m2),
-            np.full(N, m3),
+            torch.full((N,), m0),
+            torch.full((N,), m1),
+            torch.full((N,), m2),
+            torch.full((N,), m3),
         ]
     else:
         inputs13 = [
@@ -199,16 +211,18 @@ def test_phase_space_via_two_2body(mapping_and_args, rng):
             r[:, 2],
             r[:, 3],
             r[:, 4],
-            np.full(N, m0),
-            np.full(N, m1),
-            np.full(N, m2),
-            np.full(N, m3),
+            torch.full((N,), m0),
+            torch.full((N,), m1),
+            torch.full((N,), m2),
+            torch.full((N,), m3),
             p0,
         ]
     p1a, p2a, p3a, det = mapping_13.map_forward(inputs13, [])
 
     # Total volume match (mean Jacobian)
-    error1 = np.sqrt(np.var(det) / N)
-    error2 = np.sqrt(np.var(J_ref) / N)
+    error1 = math.sqrt(det.var(correction=0).item() / N)
+    error2 = math.sqrt(J_ref.var(correction=0).item() / N)
     print(f"Mean det error: {error1:.3e}, ref: {error2:.3e}")
-    assert np.mean(det) == approx(np.mean(J_ref), abs=3 * max(error1, error2))
+    assert det.mean().item() == pytest.approx(
+        J_ref.mean().item(), abs=3 * max(error1, error2)
+    )
diff --git a/madspace/tests/test_two_particle.py b/madspace/tests/test_two_particle.py
index 5d3a0eb4f..1f333d0ae 100644
--- a/madspace/tests/test_two_particle.py
+++ b/madspace/tests/test_two_particle.py
@@ -1,15 +1,21 @@
+import math
 from collections import namedtuple
 
-import numpy as np
 import pytest
-from pytest import approx
+import torch
+from madspace_testing import generator, rand, randn
 
 import madspace as ms
 
 
 @pytest.fixture
 def rng():
-    return np.random.default_rng(1234)
+    return generator(1234)
+
+
+def assert_approx(actual, expected):
+    # the default tolerances of pytest.approx
+    torch.testing.assert_close(actual, expected, rtol=1e-6, atol=1e-12)
 
 
 @pytest.fixture(
@@ -32,13 +38,13 @@ def rng():
 )
 def mapping_and_args(request):
     com = request.param["com"]
-    zeros = np.zeros(N)
+    zeros = torch.zeros(N)
     if request.param["decay"]:
         mapping = ms.TwoBodyDecay(com=com)
         if com:
 
             def make_args(point):
-                p0 = np.stack([point.m0, zeros, zeros, zeros], axis=1)
+                p0 = torch.stack([point.m0, zeros, zeros, zeros], dim=1)
                 return ([point.r1, point.r2, point.m0, point.m1, point.m2], [], p0)
 
         else:
@@ -58,9 +64,9 @@ def make_args(point):
 
             def make_args(point):
                 e0 = point.m0 / 2
-                p0 = np.stack([point.m0, zeros, zeros, zeros], axis=1)
-                pa = np.stack([e0, zeros, zeros, e0], axis=1)
-                pb = np.stack([e0, zeros, zeros, -e0], axis=1)
+                p0 = torch.stack([point.m0, zeros, zeros, zeros], dim=1)
+                pa = torch.stack([e0, zeros, zeros, e0], dim=1)
+                pb = torch.stack([e0, zeros, zeros, -e0], dim=1)
                 return [point.r1, point.r2, point.m1, point.m2], [pa, pb, *cut], p0
 
         else:
@@ -76,20 +82,24 @@ def make_args(point):
 
 
 def mass(momentum):
-    return np.sqrt(momentum[:, 0] ** 2 - np.sum(momentum[:, 1:] ** 2, axis=1))
+    return torch.sqrt(momentum[:, 0] ** 2 - torch.sum(momentum[:, 1:] ** 2, dim=1))
 
 
 InputPoint = namedtuple("InputPoint", ["r1", "r2", "m0", "m1", "m2", "p0", "pa", "pb"])
 
 N = 10_000
 
-ZEROS = np.zeros(N)
-PA = np.full((N, 4), [np.sqrt(100**2 + 200**2 + 300**2), 100.0, 200.0, 300.0])
-PB = np.full((N, 4), [np.sqrt(100**2 + 200**2 + 400**2), 100.0, 200.0, -400.0])
+ZEROS = torch.zeros(N)
+PA = torch.tensor(
+    [math.sqrt(100**2 + 200**2 + 300**2), 100.0, 200.0, 300.0]
+).repeat(N, 1)
+PB = torch.tensor(
+    [math.sqrt(100**2 + 200**2 + 400**2), 100.0, 200.0, -400.0]
+).repeat(N, 1)
 P0 = PA + PB
 M0 = mass(P0)
-M1 = np.full(N, 20.0)
-M2 = np.full(N, 30.0)
+M1 = torch.full((N,), 20.0)
+M2 = torch.full((N,), 30.0)
 
 
 @pytest.fixture(
@@ -105,32 +115,35 @@ def mass(momentum):
     ],
 )
 def fixed_input_points(request, rng):
-    r1 = rng.random(N)
-    r2 = rng.random(N)
+    r1 = rand(N, generator=rng)
+    r2 = rand(N, generator=rng)
     return InputPoint(r1, r2, *request.param)
 
 
 @pytest.fixture
 def input_points(rng):
     max_mass = 125.0
-    r1 = rng.random(N)
-    r2 = rng.random(N)
+    r1 = rand(N, generator=rng)
+    r2 = rand(N, generator=rng)
+
+    def uniform(low, high):
+        return low + (high - low) * rand(N, generator=rng)
 
-    pa = np.array([[0, 0, 4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-    ma = rng.uniform(2.5 * max_mass, 1000.0, N)
-    ea = np.sqrt(np.sum(pa**2, axis=1) + ma**2)
-    pa = np.concatenate([ea[:, None], pa], axis=1)
+    pa = torch.tensor([[0, 0, 4000.0]]) + 500.0 * randn((N, 3), generator=rng)
+    ma = uniform(2.5 * max_mass, 1000.0)
+    ea = torch.sqrt(torch.sum(pa**2, dim=1) + ma**2)
+    pa = torch.cat([ea[:, None], pa], dim=1)
 
-    pb = np.array([[0, 0, -4000.0]]) + rng.normal(0.0, 500.0, (N, 3))
-    mb = rng.uniform(20.0, 500.0, N)
-    eb = np.sqrt(np.sum(pb**2, axis=1) + mb**2)
-    pb = np.concatenate([eb[:, None], pb], axis=1)
+    pb = torch.tensor([[0, 0, -4000.0]]) + 500.0 * randn((N, 3), generator=rng)
+    mb = uniform(20.0, 500.0)
+    eb = torch.sqrt(torch.sum(pb**2, dim=1) + mb**2)
+    pb = torch.cat([eb[:, None], pb], dim=1)
 
     p0 = pa + pb
     m0 = mass(p0)
 
-    m1 = rng.uniform(1.0, max_mass, N)
-    m2 = rng.uniform(1.0, max_mass, N)
+    m1 = uniform(1.0, max_mass)
+    m2 = uniform(1.0, max_mass)
 
     return InputPoint(r1, r2, m0, m1, m2, p0, pa, pb)
 
@@ -139,7 +152,7 @@ def test_momentum_conservation(mapping_and_args, input_points):
     mapping, make_args = mapping_and_args
     inputs, conditions, p0 = make_args(input_points)
     p1, p2, det = mapping.map_forward(inputs, conditions)
-    assert p1 + p2 == approx(p0)
+    assert_approx(p1 + p2, p0)
 
 
 def test_inverse(mapping_and_args, input_points):
@@ -147,9 +160,9 @@ def test_inverse(mapping_and_args, input_points):
     inputs, conditions, p0 = make_args(input_points)
     p1, p2, det = mapping.map_forward(inputs, conditions)
     *inv_inputs, inv_det = mapping.map_inverse([p1, p2], conditions)
-    assert inv_det == approx(1 / det)
+    assert_approx(inv_det, 1 / det)
     for inp, inv_inp in zip(inputs, inv_inputs):
-        assert inp == approx(inv_inp)
+        assert_approx(inp, inv_inp)
 
 
 def test_outgoing_masses(mapping_and_args, input_points):
@@ -159,9 +172,9 @@ def test_outgoing_masses(mapping_and_args, input_points):
     m0 = mass(p1 + p2)
     m1 = mass(p1)
     m2 = mass(p2)
-    assert m0 == approx(input_points.m0)
-    assert m1 == approx(input_points.m1)
-    assert m2 == approx(input_points.m2)
+    assert_approx(m0, input_points.m0)
+    assert_approx(m1, input_points.m1)
+    assert_approx(m2, input_points.m2)
 
 
 def test_phase_space_volume(mapping_and_args, fixed_input_points):
@@ -173,7 +186,7 @@ def test_phase_space_volume(mapping_and_args, fixed_input_points):
     m1_2 = fixed_input_points.m1**2
     m2_2 = fixed_input_points.m2**2
     phase_space_vol = (
-        np.pi / (2 * s) * np.sqrt((s - m1_2 - m2_2) ** 2 - 4 * m1_2 * m2_2)
+        math.pi / (2 * s) * torch.sqrt((s - m1_2 - m2_2) ** 2 - 4 * m1_2 * m2_2)
     )
 
-    assert np.mean(det) == approx(phase_space_vol)
+    assert_approx(det.mean().expand_as(phase_space_vol), phase_space_vol)
diff --git a/madspace/tests/test_two_to_three_arcsine.py b/madspace/tests/test_two_to_three_arcsine.py
index 90b01c8c3..0ee6b8f28 100644
--- a/madspace/tests/test_two_to_three_arcsine.py
+++ b/madspace/tests/test_two_to_three_arcsine.py
@@ -15,26 +15,27 @@
 
 import math
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand, randint
 
 import madspace as ms
 
 E_BEAM = 500.0
-PA = np.array([E_BEAM, 0.0, 0.0, E_BEAM])
-PB = np.array([E_BEAM, 0.0, 0.0, -E_BEAM])
+PA = torch.tensor([E_BEAM, 0.0, 0.0, E_BEAM])
+PB = torch.tensor([E_BEAM, 0.0, 0.0, -E_BEAM])
 # r_s23 from the middle out to the edges
-EDGE = 10.0 ** -np.arange(1, 6)
-R_S23 = np.concatenate([EDGE, [0.5], 1.0 - EDGE])
+EDGE = 10.0 ** -torch.arange(1.0, 6.0)
+R_S23 = torch.cat([EDGE, torch.tensor([0.5]), 1.0 - EDGE])
 
 
 def massless(energy, theta, phi):
-    return energy * np.array(
+    return energy * torch.tensor(
         [
             1.0,
-            np.sin(theta) * np.cos(phi),
-            np.sin(theta) * np.sin(phi),
-            np.cos(theta),
+            math.sin(theta) * math.cos(phi),
+            math.sin(theta) * math.sin(phi),
+            math.cos(theta),
         ]
     )
 
@@ -45,27 +46,28 @@ def block_weights(r_s23, arcsine, s_power=0.0, m1=40.0, m2=0.0):
     )
     n = len(r_s23)
     inputs = [
-        np.zeros(n, dtype=np.int32),
-        np.asarray(r_s23, dtype=np.float64),
-        np.full(n, 0.61),
-        np.full(n, m1),
-        np.full(n, m2),
+        torch.zeros(n, dtype=torch.int32),
+        r_s23,
+        torch.full((n,), 0.61),
+        torch.full((n,), m1),
+        torch.full((n,), m2),
     ]
     conditions = [
-        np.tile(PA, (n, 1)),
-        np.tile(PB, (n, 1)),
-        np.tile(massless(350.0, 0.7, 0.3), (n, 1)),
+        PA.repeat(n, 1),
+        PB.repeat(n, 1),
+        massless(350.0, 0.7, 0.3).repeat(n, 1),
     ]
     p1, p2, det = mapping.map_forward(inputs, conditions)
-    return mapping, inputs, conditions, np.asarray(det)
+    return mapping, inputs, conditions, det
 
 
 def test_flat_weight_is_constant():
     """Flat s23 plus the arcsine map is flat in phi, which is exactly the
     measure: at fixed t1 the weight must not depend on r_s23 at all."""
     *_, det = block_weights(R_S23, arcsine=True)
-    assert np.all(np.isfinite(det)) and np.all(det > 0)
-    assert det / det[len(EDGE)] == pytest.approx(1.0, rel=1e-5)
+    assert torch.all(torch.isfinite(det)) and torch.all(det > 0)
+    ratio = det / det[len(EDGE)]
+    torch.testing.assert_close(ratio, torch.ones_like(ratio), rtol=1e-5, atol=1e-12)
 
 
 @pytest.mark.parametrize("s_power", [0.8, 1.0])
@@ -73,37 +75,40 @@ def test_weight_bounded_with_importance_sampling(s_power):
     """With 1/s23^nu importance sampling the weight is no longer constant,
     but it must stay finite and bounded out to the edges of the range."""
     *_, det = block_weights(R_S23, arcsine=True, s_power=s_power)
-    assert np.all(np.isfinite(det)) and np.all(det > 0)
-    edges = np.concatenate([det[: len(EDGE)], det[len(EDGE) + 1 :]])
+    assert torch.all(torch.isfinite(det)) and torch.all(det > 0)
+    edges = torch.cat([det[: len(EDGE)], det[len(EDGE) + 1 :]])
     # converged to a finite limit: the last two decades agree
-    assert det[len(EDGE) - 1] == pytest.approx(det[len(EDGE) - 2], rel=1e-3)
-    assert det[-1] == pytest.approx(det[-2], rel=1e-3)
-    assert np.max(edges) / np.min(edges) < 1e3
+    assert det[len(EDGE) - 1].item() == pytest.approx(
+        det[len(EDGE) - 2].item(), rel=1e-3
+    )
+    assert det[-1].item() == pytest.approx(det[-2].item(), rel=1e-3)
+    assert torch.max(edges) / torch.min(edges) < 1e3
 
 
 def test_edge_peak_without_arcsine():
     """Without the map the weight keeps its 1/|sin(phi)| = 1/(2 sqrt(u (1 - u)))
     peak, with u = r_s23 for flat sampling."""
     *_, det = block_weights(R_S23, arcsine=False)
-    ratio = det * np.sqrt(R_S23 * (1.0 - R_S23))
-    assert ratio == pytest.approx(ratio[len(EDGE)], rel=1e-5)
+    ratio = det * torch.sqrt(R_S23 * (1.0 - R_S23))
+    torch.testing.assert_close(
+        ratio, ratio[len(EDGE)].expand_as(ratio), rtol=1e-5, atol=1e-12
+    )
     assert det[len(EDGE) - 1] / det[len(EDGE)] > 100.0
 
 
 @pytest.mark.parametrize("arcsine", [True, False], ids=["arcsine", "flat"])
 @pytest.mark.parametrize("s_power", [0.0, 0.8, 1.0])
 def test_round_trip(arcsine, s_power):
-    rng = np.random.default_rng(7)
-    r_s23 = np.concatenate([rng.random(10_000), R_S23])
+    r_s23 = torch.cat([rand(10_000, generator=generator(7)), R_S23])
     mapping, inputs, conditions, det = block_weights(
         r_s23, arcsine=arcsine, s_power=s_power
     )
     p1, p2, _ = mapping.map_forward(inputs, conditions)
     *inv_inputs, det_inv = mapping.map_inverse([p1, p2], conditions)
-    assert np.asarray(inv_inputs[1]) == pytest.approx(r_s23, abs=1e-12)
-    assert np.asarray(inv_inputs[2]) == pytest.approx(inputs[2], abs=1e-8)
-    rt = np.abs(det * np.asarray(det_inv) - 1.0)
-    assert np.max(rt) < 1e-8
+    torch.testing.assert_close(inv_inputs[1], r_s23, rtol=1e-6, atol=1e-12)
+    torch.testing.assert_close(inv_inputs[2], inputs[2], rtol=1e-6, atol=1e-8)
+    rt = torch.abs(det * det_inv - 1.0)
+    assert torch.max(rt) < 1e-8
 
 
 @pytest.mark.parametrize("s_power", [0.0, 0.8, 1.0])
@@ -112,15 +117,15 @@ def test_round_trip_at_the_edges(s_power):
     1 - u in both directions, never through s23 itself, so r_s23 is recovered
     to the precision of the momenta also within 1e-12 of the edges, where
     going through s23 left an absolute error of ~1e-8."""
-    edge = 10.0 ** -np.arange(6, 13)
-    r_s23 = np.concatenate([edge, 1.0 - edge])
+    edge = 10.0 ** -torch.arange(6.0, 13.0)
+    r_s23 = torch.cat([edge, 1.0 - edge])
     mapping, inputs, conditions, det = block_weights(
         r_s23, arcsine=True, s_power=s_power
     )
     p1, p2, _ = mapping.map_forward(inputs, conditions)
     *inv_inputs, det_inv = mapping.map_inverse([p1, p2], conditions)
-    assert np.asarray(inv_inputs[1]) == pytest.approx(r_s23, rel=1e-6, abs=1e-15)
-    assert np.max(np.abs(det * np.asarray(det_inv) - 1.0)) < 1e-9
+    torch.testing.assert_close(inv_inputs[1], r_s23, rtol=1e-6, atol=1e-15)
+    assert torch.max(torch.abs(det * det_inv - 1.0)) < 1e-9
 
 
 CUTS = dict(
@@ -137,31 +142,32 @@ def sample(color_order, arcsine, seed, n, cuts=None):
     cm_energy = 1000.0
     n_out = len(color_order) - 2
     mapping = ms.ColorOrderedMapping(color_order, arcsine_s23=arcsine, **cuts)
-    rng = np.random.default_rng(seed)
-    r = rng.random((n, mapping.random_dim()))
-    d = rng.integers(0, 2, size=(n, mapping.discrete_dim())).astype(np.int32)
-    inputs = [r[:, i].copy() for i in range(r.shape[1])]
-    inputs += [d[:, j].copy() for j in range(d.shape[1])]
-    conditions = [np.full(n, cm_energy)] + [np.zeros(n)] * n_out
+    g = generator(seed)
+    r = rand((n, mapping.random_dim()), generator=g)
+    d = randint(0, 2, (n, mapping.discrete_dim()), generator=g, dtype=torch.int32)
+    inputs = [r[:, i].contiguous() for i in range(r.shape[1])]
+    inputs += [d[:, j].contiguous() for j in range(d.shape[1])]
+    conditions = [torch.full((n,), cm_energy)] + [torch.zeros(n)] * n_out
     *momenta, det = mapping.map_forward(inputs, conditions)
     # momenta of the two beams, then of the outgoing particles
-    p = np.stack([np.asarray(k) for k in momenta], axis=1)[:, 2:]
-    det = np.asarray(det) * 2.0 ** mapping.discrete_dim()
-    ok = np.all(np.isfinite(p), axis=(1, 2)) & np.isfinite(det)
+    p = torch.stack(momenta, dim=1)[:, 2:]
+    det = det * 2.0 ** mapping.discrete_dim()
+    ok = torch.isfinite(p).all(dim=2).all(dim=1) & torch.isfinite(det)
     if cuts:
-        pt = np.hypot(p[..., 1], p[..., 2])
-        ok &= np.all(pt > cuts["pt_min"][0], axis=1)
-        eta = np.arcsinh(p[..., 3] / np.maximum(pt, 1e-300))
-        phi = np.arctan2(p[..., 2], p[..., 1])
+        pt = torch.hypot(p[..., 1], p[..., 2])
+        ok &= torch.all(pt > cuts["pt_min"][0], dim=1)
+        eta = torch.asinh(p[..., 3] / torch.clamp(pt, min=1e-300))
+        phi = torch.atan2(p[..., 2], p[..., 1])
         for i in range(n_out):
             for j in range(i + 1, n_out):
                 q = p[:, i] + p[:, j]
-                m2 = q[:, 0] ** 2 - np.sum(q[:, 1:] ** 2, axis=1)
-                dphi = np.angle(np.exp(1j * (phi[:, i] - phi[:, j])))
-                dr = np.hypot(eta[:, i] - eta[:, j], dphi)
+                m2 = q[:, 0] ** 2 - torch.sum(q[:, 1:] ** 2, dim=1)
+                dphi = phi[:, i] - phi[:, j]
+                dphi = torch.atan2(torch.sin(dphi), torch.cos(dphi))
+                dr = torch.hypot(eta[:, i] - eta[:, j], dphi)
                 ok &= m2 > cuts["m_inv_min"][i][j] ** 2
                 ok &= dr > cuts["dr_min"][i][j]
-    return np.where(ok, det, 0.0)
+    return torch.where(ok, det, 0.0)
 
 
 @pytest.mark.parametrize(
@@ -173,9 +179,10 @@ def test_color_ordered_massless_volume(color_order):
     """The analytic massless 4-body volume, (pi/2)^3 s^2 / 12, with the map."""
     volume = (math.pi / 2) ** 3 * 1000.0**4 / 12
     means = [
-        sample(color_order, True, seed, 1_000_000).mean() for seed in range(200, 204)
+        sample(color_order, True, seed, 1_000_000).mean().item()
+        for seed in range(200, 204)
     ]
-    assert np.mean(means) / volume == pytest.approx(1.0, abs=0.02)
+    assert sum(means) / len(means) / volume == pytest.approx(1.0, abs=0.02)
 
 
 def test_color_ordered_cut_volume_unchanged():
@@ -183,7 +190,10 @@ def test_color_ordered_cut_volume_unchanged():
     map only covers that part. The fiducial volume must not change."""
     order = [0, 2, 3, 4, 5, 1]
     n = 1_000_000
-    on = np.concatenate([sample(order, True, s, n, CUTS) for s in range(300, 304)])
-    off = np.concatenate([sample(order, False, s, n, CUTS) for s in range(310, 314)])
-    err = math.hypot(on.std() / math.sqrt(len(on)), off.std() / math.sqrt(len(off)))
-    assert abs(on.mean() - off.mean()) < 4.0 * err
+    on = torch.cat([sample(order, True, s, n, CUTS) for s in range(300, 304)])
+    off = torch.cat([sample(order, False, s, n, CUTS) for s in range(310, 314)])
+    err = math.hypot(
+        on.std(correction=0).item() / math.sqrt(len(on)),
+        off.std(correction=0).item() / math.sqrt(len(off)),
+    )
+    assert abs(on.mean() - off.mean()).item() < 4.0 * err
diff --git a/madspace/tests/test_two_to_three_soft_limit.py b/madspace/tests/test_two_to_three_soft_limit.py
index fdf6bab6f..282e3b65e 100644
--- a/madspace/tests/test_two_to_three_soft_limit.py
+++ b/madspace/tests/test_two_to_three_soft_limit.py
@@ -19,64 +19,65 @@
 
 import math
 
-import numpy as np
 import pytest
+import torch
+from madspace_testing import generator, rand, randint
 
 import madspace as ms
 
 E_BEAM = 500.0
-PA = np.array([E_BEAM, 0.0, 0.0, E_BEAM])
-PB = np.array([E_BEAM, 0.0, 0.0, -E_BEAM])
+PA = torch.tensor([E_BEAM, 0.0, 0.0, E_BEAM])
+PB = torch.tensor([E_BEAM, 0.0, 0.0, -E_BEAM])
 # deep enough into the soft limit to reach where the old Gram determinant
 # broke down (about 1e-8)
-DELTAS = 10.0 ** -np.arange(2, 15)
+DELTAS = 10.0 ** -torch.arange(2.0, 15.0)
 
 
 def massless(energy, theta, phi):
-    return energy * np.array(
+    return energy * torch.tensor(
         [
             1.0,
-            np.sin(theta) * np.cos(phi),
-            np.sin(theta) * np.sin(phi),
-            np.cos(theta),
+            math.sin(theta) * math.cos(phi),
+            math.sin(theta) * math.sin(phi),
+            math.cos(theta),
         ]
     )
 
 
 def lsquare(p):
-    return p[..., 0] ** 2 - np.sum(p[..., 1:] ** 2, axis=-1)
+    return p[..., 0] ** 2 - torch.sum(p[..., 1:] ** 2, dim=-1)
 
 
 def scatter(p3, m1, m2):
     mapping = ms.TwoToThreeParticleScattering(0.8, 0.0, 0.0, 0.8, 0.0, 0.0)
     n = len(m1)
     inputs = [
-        np.zeros(n, dtype=np.int32),
-        np.full(n, 0.37),
-        np.full(n, 0.61),
+        torch.zeros(n, dtype=torch.int32),
+        torch.full((n,), 0.37),
+        torch.full((n,), 0.61),
         m1,
         m2,
     ]
-    conditions = [np.tile(PA, (n, 1)), np.tile(PB, (n, 1)), p3]
+    conditions = [PA.repeat(n, 1), PB.repeat(n, 1), p3]
     _k1, _k2, det = mapping.map_forward(inputs, conditions)
-    return np.asarray(det)
+    return det
 
 
 def assert_smooth(values):
     """Every value is finite and within a factor 2 of the first one. (The s23
     sampling moves the limit by a few percent once s23 drops below the 1e-2
     GeV^2 offset of its power map, so this is not a tighter check.)"""
-    assert np.all(np.isfinite(values)), values
+    assert torch.all(torch.isfinite(values)), values
     ratio = values / values[0]
-    assert np.all((ratio > 0.5) & (ratio < 2.0)), ratio
+    assert torch.all((ratio > 0.5) & (ratio < 2.0)), ratio
 
 
 def test_jacobian_with_soft_peeled_particle():
     """m1 -> sqrt(s12) leaves particle 2 soft. Its phase space closes
     linearly, so det / delta must tend to a constant."""
-    p3 = np.tile(massless(350.0, 0.7, 0.3), (len(DELTAS), 1))
-    sqrt_s12 = math.sqrt(lsquare(PA + PB - p3[0]))
-    det = scatter(p3, sqrt_s12 * (1.0 - DELTAS), np.zeros(len(DELTAS)))
+    p3 = massless(350.0, 0.7, 0.3).repeat(len(DELTAS), 1)
+    sqrt_s12 = math.sqrt(lsquare(PA + PB - p3[0]).item())
+    det = scatter(p3, sqrt_s12 * (1.0 - DELTAS), torch.zeros(len(DELTAS)))
     assert_smooth(det / DELTAS)
 
 
@@ -84,9 +85,9 @@ def test_jacobian_with_soft_third_particle():
     """A soft third momentum (the particle peeled in the step before) leaves
     the rest of the kinematics regular, so the Jacobian must tend to a
     constant."""
-    p3 = np.stack([massless(350.0 * delta, 0.7, 0.3) for delta in DELTAS])
-    sqrt_s12 = np.sqrt(lsquare(PA + PB - p3))
-    det = scatter(p3, 0.4 * sqrt_s12, np.zeros(len(DELTAS)))
+    p3 = torch.stack([massless(350.0 * delta, 0.7, 0.3) for delta in DELTAS.tolist()])
+    sqrt_s12 = torch.sqrt(lsquare(PA + PB - p3))
+    det = scatter(p3, 0.4 * sqrt_s12, torch.zeros(len(DELTAS)))
     assert_smooth(det)
 
 
@@ -110,21 +111,21 @@ def test_color_ordered_massless_volume(color_order):
     n = 2_000_000
     means = []
     for seed in range(100, 108):
-        rng = np.random.default_rng(seed)
-        r = rng.random((n, mapping.random_dim()))
-        d = rng.integers(0, 2, size=(n, mapping.discrete_dim())).astype(np.int32)
-        inputs = [r[:, i].copy() for i in range(r.shape[1])]
-        inputs += [d[:, j].copy() for j in range(d.shape[1])]
-        conditions = [np.full(n, cm_energy)] + [np.zeros(n)] * n_out
+        g = generator(seed)
+        r = rand((n, mapping.random_dim()), generator=g)
+        d = randint(0, 2, (n, mapping.discrete_dim()), generator=g, dtype=torch.int32)
+        inputs = [r[:, i].contiguous() for i in range(r.shape[1])]
+        inputs += [d[:, j].contiguous() for j in range(d.shape[1])]
+        conditions = [torch.full((n,), cm_energy)] + [torch.zeros(n)] * n_out
         *momenta, det = mapping.map_forward(inputs, conditions)
-        p = np.stack([np.asarray(k) for k in momenta], axis=1)
-        det = np.asarray(det) * 2.0 ** mapping.discrete_dim()
-        det = np.where(np.all(np.isfinite(p), axis=(1, 2)), det, 0.0)
-        means.append(np.nan_to_num(det).mean())
+        p = torch.stack(momenta, dim=1)
+        det = det * 2.0 ** mapping.discrete_dim()
+        det = torch.where(torch.isfinite(p).all(dim=2).all(dim=1), det, 0.0)
+        means.append(torch.nan_to_num(det).mean().item())
     # The weights keep an integrable 1/sin(phi) peak at the edges of the s23
     # range, so the error estimate itself is noisy; a fixed tolerance well
     # above the spread (about 1%) is the robust check.
-    assert np.mean(means) / volume == pytest.approx(1.0, abs=0.05)
+    assert sum(means) / len(means) / volume == pytest.approx(1.0, abs=0.05)
 
 
 def extreme_points():
@@ -147,7 +148,8 @@ def extreme_points():
             for r_s23 in corners:
                 for r_t1 in corners:
                     for index in (0, 1):
-                        points.append((p3, m1(lsquare(PA + PB - p3)), r_s23, r_t1, index))
+                        s12 = lsquare(PA + PB - p3).item()
+                        points.append((p3, m1(s12), r_s23, r_t1, index))
     return points
 
 
@@ -165,17 +167,20 @@ def test_finite_at_extreme_points(arcsine, s_power):
         0.8, 0.0, 0.0, s_power, 0.0, 0.0, arcsine_s23=arcsine
     )
     inputs = [
-        np.array([pt[4] for pt in points], dtype=np.int32),
-        np.array([pt[2] for pt in points]),
-        np.array([pt[3] for pt in points]),
-        np.array([pt[1] for pt in points]),
-        np.zeros(n),
+        torch.tensor([pt[4] for pt in points], dtype=torch.int32),
+        torch.tensor([pt[2] for pt in points]),
+        torch.tensor([pt[3] for pt in points]),
+        torch.tensor([pt[1] for pt in points]),
+        torch.zeros(n),
+    ]
+    conditions = [
+        PA.repeat(n, 1),
+        PB.repeat(n, 1),
+        torch.stack([pt[0] for pt in points]),
     ]
-    conditions = [np.tile(PA, (n, 1)), np.tile(PB, (n, 1)), np.stack([pt[0] for pt in points])]
     p1, p2, det = mapping.map_forward(inputs, conditions)
-    p1, p2, det = np.asarray(p1), np.asarray(p2), np.asarray(det)
-    assert np.all(np.isfinite(p1)) and np.all(np.isfinite(p2))
-    assert np.all(np.isfinite(det)) and np.all(det >= 0)
+    assert torch.all(torch.isfinite(p1)) and torch.all(torch.isfinite(p2))
+    assert torch.all(torch.isfinite(det)) and torch.all(det >= 0)
     # On shell where the weight is not zero: the soft-p3 points used to be off
     # by ~1e4 GeV^2, |p^2 - m^2| / E^2 ~ 1e-2. What is left comes from s12 -> 0,
     # a pair collinear to ~1e-6 whose mass pa + pb - p3 fixes only to ~1e-4,
@@ -183,11 +188,11 @@ def test_finite_at_extreme_points(arcsine, s_power):
     # sqrt(eps) of its energy, which momenta cannot carry at all.)
     live = det > 0
     m1 = inputs[3][live]
-    off1 = np.abs(lsquare(p1[live]) - m1**2) / p1[live, 0] ** 2
-    off2 = np.abs(lsquare(p2[live])) / p2[live, 0] ** 2
-    assert np.max(off1) < 1e-8 and np.max(off2) < 1e-8
+    off1 = torch.abs(lsquare(p1[live]) - m1**2) / p1[live, 0] ** 2
+    off2 = torch.abs(lsquare(p2[live])) / p2[live, 0] ** 2
+    assert torch.max(off1) < 1e-8 and torch.max(off2) < 1e-8
     for out in mapping.map_inverse([p1, p2], conditions):
-        assert np.all(np.isfinite(np.asarray(out)))
+        assert torch.all(torch.isfinite(out))
 
 
 def test_arcsine_rejects_breit_wigner():