diff --git a/benchmarks/ubp.py b/benchmarks/ubp.py
new file mode 100644
index 00000000..8e16dfb5
--- /dev/null
+++ b/benchmarks/ubp.py
@@ -0,0 +1,235 @@
+import math
+from typing import cast
+
+import diffrax
+import equinox as eqx
+import equinox.internal as eqxi
+import jax
+import jax.numpy as jnp
+import jax.random as jr
+import jax.tree_util as jtu
+import lineax.internal as lxi
+from jaxtyping import PRNGKeyArray, PyTree
+from lineax.internal import complex_to_real_dtype
+
+
+class OldBrownianPath(diffrax.AbstractBrownianPath):
+    shape: PyTree[jax.ShapeDtypeStruct] = eqx.field(static=True)
+    levy_area: type[
+        diffrax.BrownianIncrement
+        | diffrax.SpaceTimeLevyArea
+        | diffrax.SpaceTimeTimeLevyArea
+    ] = eqx.field(static=True)
+    key: PRNGKeyArray
+
+    def __init__(
+        self,
+        shape,
+        key,
+        levy_area=diffrax.BrownianIncrement,
+    ):
+        self.shape = (
+            jax.ShapeDtypeStruct(shape, lxi.default_floating_dtype())
+            if diffrax._misc.is_tuple_of_ints(shape)
+            else shape
+        )
+        self.key = key
+        self.levy_area = levy_area
+
+        if any(
+            not jnp.issubdtype(x.dtype, jnp.inexact)
+            for x in jtu.tree_leaves(self.shape)
+        ):
+            raise ValueError("OldBrownianPath dtypes all have to be floating-point.")
+
+    @property
+    def t0(self):
+        return -jnp.inf
+
+    @property
+    def t1(self):
+        return jnp.inf
+
+
+    def __call__(
+        self,
+        t0,
+        brownian_state,
+        t1=None,
+        left=True,
+        use_levy=False,
+    ):
+        return self.evaluate(t0, t1, left, use_levy), brownian_state
+
+    @eqx.filter_jit
+    def evaluate(
+        self,
+        t0,
+        t1=None,
+        left=True,
+        use_levy=False,
+        index=None
+    ):
+        del left, index
+        if t1 is None:
+            dtype = jnp.result_type(t0)
+            t1 = t0
+            t0 = jnp.array(0, dtype)
+        else:
+            with jax.numpy_dtype_promotion("standard"):
+                dtype = jnp.result_type(t0, t1)
+            t0 = jnp.astype(t0, dtype)
+            t1 = jnp.astype(t1, dtype)
+        t0 = eqxi.nondifferentiable(t0, name="t0")
+        t1 = eqxi.nondifferentiable(t1, name="t1")
+        t1 = cast(diffrax._custom_types.RealScalarLike, t1)
+        t0_ = diffrax._misc.force_bitcast_convert_type(t0, jnp.int32)
+        t1_ = diffrax._misc.force_bitcast_convert_type(t1, jnp.int32)
+        key = jr.fold_in(self.key, t0_)
+        key = jr.fold_in(key, t1_)
+        key = diffrax._misc.split_by_tree(key, self.shape)
+        out = jtu.tree_map(
+            lambda key, shape: self._evaluate_leaf(
+                t0, t1, key, shape, self.levy_area, use_levy
+            ),
+            key,
+            self.shape,
+        )
+        if use_levy:
+            out = diffrax._custom_types.levy_tree_transpose(self.shape, out)
+            assert isinstance(out, self.levy_area)
+        return out
+
+    @staticmethod
+    def _evaluate_leaf(
+        t0,
+        t1,
+        key,
+        shape,
+        levy_area,
+        use_levy,
+    ):
+        w_std = jnp.sqrt(t1 - t0).astype(shape.dtype)
+        dt = jnp.asarray(t1 - t0, dtype=complex_to_real_dtype(shape.dtype))
+
+        if levy_area is diffrax.SpaceTimeTimeLevyArea:
+            key_w, key_hh, key_kk = jr.split(key, 3)
+            w = jr.normal(key_w, shape.shape, shape.dtype) * w_std
+            hh_std = w_std / math.sqrt(12)
+            hh = jr.normal(key_hh, shape.shape, shape.dtype) * hh_std
+            kk_std = w_std / math.sqrt(720)
+            kk = jr.normal(key_kk, shape.shape, shape.dtype) * kk_std
+            levy_val = diffrax.SpaceTimeTimeLevyArea(dt=dt, W=w, H=hh, K=kk)
+
+        elif levy_area is diffrax.SpaceTimeLevyArea:
+            key_w, key_hh = jr.split(key, 2)
+            w = jr.normal(key_w, shape.shape, shape.dtype) * w_std
+            hh_std = w_std / math.sqrt(12)
+            hh = jr.normal(key_hh, shape.shape, shape.dtype) * hh_std
+            levy_val = diffrax.SpaceTimeLevyArea(dt=dt, W=w, H=hh)
+        elif levy_area is diffrax.BrownianIncrement:
+            w = jr.normal(key, shape.shape, shape.dtype) * w_std
+            levy_val = diffrax.BrownianIncrement(dt=dt, W=w)
+        else:
+            assert False
+
+        if use_levy:
+            return levy_val
+        return w
+
+
+# https://github.com/patrick-kidger/diffrax/issues/517
+key = jax.random.key(42)
+# t0 = 0
+# t1 = 100
+# y0 = 1.0
+# ndt = 4000
+# dt = (t1 - t0) / (ndt - 1)
+# drift = lambda t, y, args: -y
+# diffusion = lambda t, y, args: 0.2
+t0 = 0
+t1 = 1
+y0 = 1.0
+ndt = 40010
+dt = (t1 - t0) / (ndt - 1)
+drift = lambda t, y, args: -y
+diffusion = lambda t, y, args: 0.2
+# saveat = diffrax.SaveAt(ts=jnp.linspace(t0, t1, ndt))
+saveat = diffrax.SaveAt(steps=True)
+
+brownian_motion = diffrax.VirtualBrownianTree(t0, t1, tol=1e-3, shape=(), key=key)
+ubp = OldBrownianPath(shape=(), key=key)
+new_ubp_pre = diffrax.UnsafeBrownianPath(shape=(), key=key, precompute=ndt + 10)
+
+solver = diffrax.Euler()
+
+terms = diffrax.MultiTerm(
+    diffrax.ODETerm(drift), diffrax.ControlTerm(diffusion, brownian_motion)
+)
+terms_old = diffrax.MultiTerm(
+    diffrax.ODETerm(drift), diffrax.ControlTerm(diffusion, ubp)
+)
+terms_new_precompute = diffrax.MultiTerm(
+    diffrax.ODETerm(drift), diffrax.ControlTerm(diffusion, new_ubp_pre)
+)
+
+
+@jax.jit
+def diffrax_vbt():
+    return diffrax.diffeqsolve(
+        terms, solver, t0, t1, dt0=dt, y0=y0, saveat=saveat, throw=False
+    ).ys
+
+
+@jax.jit
+def diffrax_old():
+    return diffrax.diffeqsolve(
+        terms_old, solver, t0, t1, dt0=dt, y0=y0, saveat=saveat, throw=False
+    ).ys
+
+
+
+@jax.jit
+def diffrax_new_pre():
+    return diffrax.diffeqsolve(
+        terms_new_precompute, solver, t0, t1, dt0=dt, y0=y0, saveat=saveat, throw=False
+    ).ys
+
+
+@jax.jit
+def homemade_simu():
+    dWs = jnp.sqrt(dt) * jax.random.normal(key, (ndt,))
+
+    def step(y, dW):
+        dy = drift(None, y, None) * dt + diffusion(None, y, None) * dW
+        return y + dy, y
+
+    return jax.lax.scan(step, y0, dWs)[-1]
+
+
+_ = diffrax_vbt().block_until_ready()
+_ = diffrax_old().block_until_ready()
+_ = diffrax_new_pre().block_until_ready()
+_ = homemade_simu().block_until_ready()
+
+from timeit import Timer
+
+
+num_runs = 10
+
+timer = Timer(stmt="_ = diffrax_vbt().block_until_ready()", globals=globals())
+total_time = timer.timeit(number=num_runs)
+print(f"VBT: {total_time / num_runs:.6f}")
+
+timer = Timer(stmt="_ = diffrax_old().block_until_ready()", globals=globals())
+total_time = timer.timeit(number=num_runs)
+print(f"Old UBP: {total_time / num_runs:.6f}")
+
+
+timer = Timer(stmt="_ = diffrax_new_pre().block_until_ready()", globals=globals())
+total_time = timer.timeit(number=num_runs)
+print(f"New UBP + Precompute: {total_time / num_runs:.6f}")
+
+timer = Timer(stmt="_ = homemade_simu().block_until_ready()", globals=globals())
+total_time = timer.timeit(number=num_runs)
+print(f"Pure Jax: {total_time / num_runs:.6f}")
\ No newline at end of file
diff --git a/diffrax/_adjoint.py b/diffrax/_adjoint.py
index 7bc081b9..a1e5c99e 100644
--- a/diffrax/_adjoint.py
+++ b/diffrax/_adjoint.py
@@ -274,12 +274,12 @@ def loop(
         **kwargs,
     ):
         del throw, passed_solver_state, passed_controller_state
-        if is_unsafe_sde(terms):
-            raise ValueError(
-                "`adjoint=RecursiveCheckpointAdjoint()` does not support "
-                "`UnsafeBrownianPath`. Consider using `adjoint=ForwardMode()` "
-                "instead."
-            )
+        # if is_unsafe_sde(terms):
+        #     raise ValueError(
+        #         "`adjoint=RecursiveCheckpointAdjoint()` does not support "
+        #         "`UnsafeBrownianPath`. Consider using `adjoint=ForwardMode()` "
+        #         "instead."
+        #     )
         if self.checkpoints is None and max_steps is None:
             inner_while_loop = ft.partial(_inner_loop, kind="lax")
             outer_while_loop = ft.partial(_outer_loop, kind="lax")
diff --git a/diffrax/_brownian/base.py b/diffrax/_brownian/base.py
index e18ff623..9ca2e232 100644
--- a/diffrax/_brownian/base.py
+++ b/diffrax/_brownian/base.py
@@ -11,7 +11,7 @@
     SpaceTimeLevyArea,
 )
 from .._path import AbstractPath
-
+from jaxtyping import Int, Array
 
 _Control = TypeVar("_Control", bound=PyTree[Array] | AbstractBrownianIncrement)
 
@@ -28,6 +28,7 @@ def evaluate(
         t1: RealScalarLike | None = None,
         left: bool = True,
         use_levy: bool = False,
+        index: None | Int[Array, ""] = None
     ) -> _Control:
         r"""Samples a Brownian increment $w(t_1) - w(t_0)$.
 
diff --git a/diffrax/_brownian/path.py b/diffrax/_brownian/path.py
index 61f49644..650baf8a 100644
--- a/diffrax/_brownian/path.py
+++ b/diffrax/_brownian/path.py
@@ -8,7 +8,7 @@
 import jax.random as jr
 import jax.tree_util as jtu
 import lineax.internal as lxi
-from jaxtyping import Array, PRNGKeyArray, PyTree
+from jaxtyping import Array, PRNGKeyArray, PyTree, Int, Float
 from lineax.internal import complex_to_real_dtype
 
 from .._custom_types import (
@@ -26,7 +26,6 @@
 )
 from .base import AbstractBrownianPath
 
-
 class UnsafeBrownianPath(AbstractBrownianPath):
     """Brownian simulation that is only suitable for certain cases.
 
@@ -66,6 +65,7 @@ class UnsafeBrownianPath(AbstractBrownianPath):
         eqx.field(static=True)
     )
     key: PRNGKeyArray
+    arr: PyTree[Float[Array, " steps"]] | None
 
     def __init__(
         self,
@@ -74,6 +74,7 @@ def __init__(
         levy_area: type[
             BrownianIncrement | SpaceTimeLevyArea | SpaceTimeTimeLevyArea
         ] = BrownianIncrement,
+        precompute: int | None = None
     ):
         """**Arguments:**
 
@@ -91,8 +92,19 @@ def __init__(
             if is_tuple_of_ints(shape)
             else shape
         )
+        key, subkey = jax.random.split(key)
         self.key = key
         self.levy_area = levy_area
+        # seems bad to define precompute and max_steps, should be the same I imagine? make eqx.tree_at in integrate? Idk
+        if precompute is None:
+            self.arr = None
+        else:
+            subkeys = split_by_tree(subkey, self.shape)
+            self.arr = jax.tree.map(
+                lambda subkey, shape: self._generate_noise(subkey, shape, precompute),
+                subkeys,
+                self.shape,
+            )
 
         if any(
             not jnp.issubdtype(x.dtype, jnp.inexact)
@@ -100,6 +112,23 @@ def __init__(
         ):
             raise ValueError("UnsafeBrownianPath dtypes all have to be floating-point.")
 
+    def _generate_noise(
+        self,
+        key: PRNGKeyArray,
+        shape: jax.ShapeDtypeStruct,
+        max_steps: int,
+    ) -> Float[Array, "..."]:
+        if self.levy_area is SpaceTimeTimeLevyArea:
+            noise = jr.normal(key, (max_steps, 3, *shape.shape), shape.dtype)
+        elif self.levy_area is SpaceTimeLevyArea:
+            noise = jr.normal(key, (max_steps, 2, *shape.shape), shape.dtype)
+        elif self.levy_area is BrownianIncrement:
+            noise = jr.normal(key, (max_steps, *shape.shape), shape.dtype)
+        else:
+            assert False
+
+        return noise
+
     @property
     def t0(self):
         return -jnp.inf
@@ -115,6 +144,7 @@ def evaluate(
         t1: RealScalarLike | None = None,
         left: bool = True,
         use_levy: bool = False,
+        index: None | Int[Array, ""] = None
     ) -> PyTree[Array] | AbstractBrownianIncrement:
         """Implements [`diffrax.AbstractBrownianPath.evaluate`][]."""
         del left
@@ -127,6 +157,20 @@ def evaluate(
                 dtype = jnp.result_type(t0, t1)
             t0 = jnp.astype(t0, dtype)
             t1 = jnp.astype(t1, dtype)
+
+        if self.arr is not None:
+            out = jax.tree.map(
+                lambda shape, noise: self._evaluate_leaf_precomputed(
+                    t0, t1, shape, self.levy_area, use_levy, noise
+                ),
+                self.shape,
+                jax.tree.map(lambda x: x[index], self.arr),
+            )
+            if use_levy:
+                out = levy_tree_transpose(self.shape, out)
+                assert isinstance(out, self.levy_area)
+            return out
+
         t0 = eqxi.nondifferentiable(t0, name="t0")
         t1 = eqxi.nondifferentiable(t1, name="t1")
         t1 = cast(RealScalarLike, t1)
@@ -183,3 +227,38 @@ def _evaluate_leaf(
         if use_levy:
             return levy_val
         return w
+
+    @staticmethod
+    def _evaluate_leaf_precomputed(
+        t0: RealScalarLike,
+        t1: RealScalarLike,
+        shape: jax.ShapeDtypeStruct,
+        levy_area: type[BrownianIncrement | SpaceTimeLevyArea | SpaceTimeTimeLevyArea],
+        use_levy: bool,
+        noises: Float[Array, "..."],
+    ):
+        w_std = jnp.sqrt(t1 - t0).astype(shape.dtype)
+        dt = jnp.asarray(t1 - t0, dtype=complex_to_real_dtype(shape.dtype))
+
+        if levy_area is SpaceTimeTimeLevyArea:
+            w = noises[0] * w_std
+            hh_std = w_std / math.sqrt(12)
+            hh = noises[1] * hh_std
+            kk_std = w_std / math.sqrt(720)
+            kk = noises[2] * kk_std
+            levy_val = SpaceTimeTimeLevyArea(dt=dt, W=w, H=hh, K=kk)
+
+        elif levy_area is SpaceTimeLevyArea:
+            w = noises[0] * w_std
+            hh_std = w_std / math.sqrt(12)
+            hh = noises[1] * hh_std
+            levy_val = SpaceTimeLevyArea(dt=dt, W=w, H=hh)
+        elif levy_area is BrownianIncrement:
+            w = noises * w_std
+            levy_val = BrownianIncrement(dt=dt, W=w)
+        else:
+            assert False
+
+        if use_levy:
+            return levy_val
+        return w
\ No newline at end of file
diff --git a/diffrax/_brownian/tree.py b/diffrax/_brownian/tree.py
index 8a430668..9001e0c8 100644
--- a/diffrax/_brownian/tree.py
+++ b/diffrax/_brownian/tree.py
@@ -9,7 +9,7 @@
 import jax.random as jr
 import jax.tree_util as jtu
 import lineax.internal as lxi
-from jaxtyping import Array, Inexact, PRNGKeyArray, PyTree
+from jaxtyping import Array, Inexact, PRNGKeyArray, PyTree, Int
 from lineax.internal import complex_to_real_dtype
 
 from .._custom_types import (
@@ -329,6 +329,7 @@ def evaluate(
         t1: RealScalarLike | None = None,
         left: bool = True,
         use_levy: bool = False,
+        index: None | Int[Array, ""] = None
     ) -> PyTree[Array] | AbstractBrownianIncrement:
         """Implements [`diffrax.AbstractBrownianPath.evaluate`][]."""
         del left
diff --git a/diffrax/_integrate.py b/diffrax/_integrate.py
index 5441e0a9..4700ed90 100644
--- a/diffrax/_integrate.py
+++ b/diffrax/_integrate.py
@@ -378,6 +378,7 @@ def body_fun_aux(state):
             args,
             state.solver_state,
             state.made_jump,
+            state.num_steps
         )
 
         # e.g. if someone has a sqrt(y) in the vector field, and dt0 is so large that
diff --git a/diffrax/_solver/base.py b/diffrax/_solver/base.py
index d4a476e3..ae28394b 100644
--- a/diffrax/_solver/base.py
+++ b/diffrax/_solver/base.py
@@ -22,7 +22,7 @@
 else:
     from equinox import AbstractClassVar, AbstractVar
 from equinox.internal import ω
-from jaxtyping import PyTree
+from jaxtyping import PyTree, Int, Array
 
 from .._custom_types import Args, BoolScalarLike, DenseInfo, RealScalarLike, VF, Y
 from .._heuristics import is_sde
@@ -147,6 +147,7 @@ def step(
         args: Args,
         solver_state: _SolverState,
         made_jump: BoolScalarLike,
+        index: Int[Array, ""] | None
     ) -> tuple[Y, Y | None, DenseInfo, _SolverState, RESULTS]:
         """Make a single step of the solver.
 
diff --git a/diffrax/_solver/euler.py b/diffrax/_solver/euler.py
index b1a323f7..968adfaf 100644
--- a/diffrax/_solver/euler.py
+++ b/diffrax/_solver/euler.py
@@ -8,7 +8,7 @@
 from .._solution import RESULTS
 from .._term import AbstractTerm
 from .base import AbstractItoSolver
-
+from jaxtyping import Int, Array
 
 _ErrorEstimate: TypeAlias = None
 _SolverState: TypeAlias = None
@@ -53,9 +53,10 @@ def step(
         args: Args,
         solver_state: _SolverState,
         made_jump: BoolScalarLike,
+        index: Int[Array, ""] | None
     ) -> tuple[Y, _ErrorEstimate, DenseInfo, _SolverState, RESULTS]:
         del solver_state, made_jump
-        control = terms.contr(t0, t1)
+        control = terms.contr(t0, t1, index)
         y1 = (y0**ω + terms.vf_prod(t0, y0, args, control) ** ω).ω
         dense_info = dict(y0=y0, y1=y1)
         return y1, None, dense_info, None, RESULTS.successful
diff --git a/diffrax/_term.py b/diffrax/_term.py
index ab231b2e..536cbd85 100644
--- a/diffrax/_term.py
+++ b/diffrax/_term.py
@@ -62,7 +62,13 @@ def vf(self, t: RealScalarLike, y: Y, args: Args) -> _VF:
         pass
 
     @abc.abstractmethod
-    def contr(self, t0: RealScalarLike, t1: RealScalarLike, **kwargs) -> _Control:
+    def contr(
+        self,
+        t0: RealScalarLike,
+        t1: RealScalarLike,
+        index: IntScalarLike | None = None,
+        **kwargs,
+    ) -> _Control:
         r"""The control.
 
         Represents the $\mathrm{d}t$ in an ODE, or the $\mathrm{d}w(t)$ in an SDE, etc.
@@ -210,7 +216,13 @@ def _broadcast_and_upcast(oi, yi):
 
         return jtu.tree_map(_broadcast_and_upcast, out, y)
 
-    def contr(self, t0: RealScalarLike, t1: RealScalarLike, **kwargs) -> RealScalarLike:
+    def contr(
+        self,
+        t0: RealScalarLike,
+        t1: RealScalarLike,
+        index: IntScalarLike | None = None,
+        **kwargs,
+    ) -> RealScalarLike:
         return t1 - t0
 
     def prod(self, vf: _VF, control: RealScalarLike) -> Y:
@@ -417,8 +429,14 @@ def __init__(
     def vf(self, t: RealScalarLike, y: Y, args: Args) -> VF:
         return self.vector_field(t, y, args)
 
-    def contr(self, t0: RealScalarLike, t1: RealScalarLike, **kwargs) -> _Control:
-        return self.control.evaluate(t0, t1, **kwargs)
+    def contr(
+        self,
+        t0: RealScalarLike,
+        t1: RealScalarLike,
+        index: IntScalarLike | None = None,
+        **kwargs,
+    ) -> _Control:
+        return self.control.evaluate(t0, t1, index=index, **kwargs)
 
     def prod(self, vf: _VF, control: _Control) -> Y:
         if isinstance(vf, lx.AbstractLinearOperator):
@@ -695,9 +713,15 @@ def vf(self, t: RealScalarLike, y: Y, args: Args) -> tuple[PyTree[ArrayLike], ..
         return tuple(term.vf(t, y, args) for term in self.terms)
 
     def contr(
-        self, t0: RealScalarLike, t1: RealScalarLike, **kwargs
+        self,
+        t0: RealScalarLike,
+        t1: RealScalarLike,
+        index: IntScalarLike | None = None,
+        **kwargs,
     ) -> tuple[PyTree[ArrayLike], ...]:
-        return tuple(term.contr(t0, t1, **kwargs) for term in self.terms)
+        return tuple(
+            term.contr(t0, t1, index, **kwargs) for term in self.terms
+        )
 
     def prod(
         self, vf: tuple[PyTree[ArrayLike], ...], control: tuple[PyTree[ArrayLike], ...]
@@ -740,10 +764,10 @@ def vf(self, t: RealScalarLike, y: Y, args: Args) -> _VF:
             t = t * self.direction
         return self.term.vf(t, y, args)
 
-    def contr(self, t0: RealScalarLike, t1: RealScalarLike, **kwargs) -> _Control:
+    def contr(self, t0: RealScalarLike, t1: RealScalarLike, index: IntScalarLike | None = None, **kwargs) -> _Control:
         _t0 = jnp.where(self.direction == 1, t0, -t1)
         _t1 = jnp.where(self.direction == 1, t1, -t0)
-        return (self.direction * self.term.contr(_t0, _t1, **kwargs) ** ω).ω
+        return (self.direction * self.term.contr(_t0, _t1, index, **kwargs) ** ω).ω
 
     def prod(self, vf: _VF, control: _Control) -> Y:
         with jax.numpy_dtype_promotion("standard"):