How to use reversible integrators?

I wanted a time-reversible integrator, and according to the documentation SemiImplicitEuler() should work.
However, I ran the following:

def f(t, p, args):
    return -p

def g(t, q, args):
    return -q

term = (diffrax.ODETerm(f), diffrax.ODETerm(g))
solver = diffrax.SemiImplicitEuler()

q0 = jnp.array([12., 3.])
p0 = jnp.array([1., 2.])

# Forward pass
y0 = (q0, p0)
solution = diffrax.diffeqsolve(term, solver, t0=0, t1=1, dt0=0.01, y0=y0,
                               saveat=diffrax.SaveAt(t1=True, solver_state=True))
y1 = jax.tree_map(lambda arr: arr[0], solution.ys)

# Backward pass
solution = diffrax.diffeqsolve(term, solver, t0=1, t1=0, dt0=-0.01, y0=y1,
                               solver_state=solution.solver_state)
y0_computed = jax.tree_map(lambda arr: arr[0], solution.ys)

I see that y0_computed is:

(Array([11.797033 ,  2.9734273]),
 Array([0.85248  , 1.9949753])

which is close to, but different from y0. Am I using the wrong arguments?

Start by running the Python reproduction in the issue with SemiImplicitEuler(), SaveAt(t1=True, solver_state=True), and the forward and backward diffeqsolve calls. Compare the returned y0_computed with y0 and determine whether the arguments are incorrect or the reversible-integrator behavior is inconsistent; the issue is done when the expected behavior is established and documented or corrected.