def fast_lyapunov_max(w0, hamiltonian, dt, n_steps, d0=1e-5,
                      n_steps_per_pullback=10, noffset_orbits=2, t1=0.,
                      atol=1E-10, rtol=1E-10, nmax=0, return_orbit=True):
    from .lyapunov import dop853_lyapunov_max, dop853_lyapunov_max_dont_save
    if isinstance(hamiltonian, PotentialBase):
        from ..potential import Hamiltonian
        hamiltonian = Hamiltonian(hamiltonian)
    if not hamiltonian.c_enabled:
        raise TypeError("Input Hamiltonian must contain a C-implemented "
                        "potential and frame.")
    if not isinstance(w0, PhaseSpacePosition):
        w0 = np.asarray(w0)
        ndim = w0.shape[0]//2
        w0 = PhaseSpacePosition(pos=w0[:ndim],
                                vel=w0[ndim:])
    _w0 = np.squeeze(w0.w(hamiltonian.units))
    if _w0.ndim > 1:
        raise ValueError("Can only compute fast Lyapunov exponent for a single orbit.")
    if return_orbit:
        t,w,l = dop853_lyapunov_max(hamiltonian, _w0,
                                    dt, n_steps+1, t1,
                                    d0, n_steps_per_pullback, noffset_orbits,
                                    atol, rtol, nmax)
        w = np.rollaxis(w, -1)
        try:
            tunit = hamiltonian.units['time']
        except (TypeError, AttributeError):
            tunit = u.dimensionless_unscaled
        orbit = Orbit.from_w(w=w, units=hamiltonian.units,
                             t=t*tunit, hamiltonian=hamiltonian)
        return l/tunit, orbit
    else:
        l = dop853_lyapunov_max_dont_save(hamiltonian, _w0,
                                          dt, n_steps+1, t1,
                                          d0, n_steps_per_pullback, noffset_orbits,
                                          atol, rtol, nmax)
        try:
            tunit = hamiltonian.units['time']
        except (TypeError, AttributeError):
            tunit = u.dimensionless_unscaled
        return l/tunit