<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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):
 """ Compute the maximum Lyapunov exponent using a C-implemented estimator that uses the DOPRI853 integrator. Parameters w0 : `~gala.dynamics.PhaseSpacePosition`, array_like Initial conditions. hamiltonian : `~gala.potential.Hamiltonian` dt : numeric Timestep. n_steps : int Number of steps to run for. d0 : numeric (optional) The initial separation. n_steps_per_pullback : int (optional) Number of steps to run before re-normalizing the offset vectors. noffset_orbits : int (optional) Number of offset orbits to run. t1 : numeric (optional) Time of initial conditions. Assumed to be t=0. return_orbit : bool (optional) Store the full orbit for the parent and all offset orbits. Returns ------- LEs : :class:`~astropy.units.Quantity` Lyapunov exponents calculated from each offset / deviation orbit. orbit : `~gala.dynamics.Orbit` (optional) """


    from .lyapunov import dop853_lyapunov_max, dop853_lyapunov_max_dont_save

    # TODO: remove in v1.0
    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