  A unified fast time-stepping method for both fractional integral and
derivative operators is proposed. The fractional operator is decomposed into a
local part with memory length $\Delta T$ and a history part, where the local
part is approximated by the direct convolution method and the history part is
approximated by a fast memory-saving method. The fast method has
$O(n_0+\sum_{\ell}^L{q}_{\alpha}(N_{\ell}))$ active memory and $O(n_0n_T+
(n_T-n_0)\sum_{\ell}^L{q}_{\alpha}(N_{\ell}))$ operations, where
$L=\log(n_T-n_0)$, $n_0={\Delta T}/\tau,n_T=T/\tau$, $\tau$ is the stepsize,
$T$ is the final time, and ${q}_{\alpha}{(N_{\ell})}$ is the number of
quadrature points used in the truncated Laguerre--Gauss (LG) quadrature. The
error bound of the present fast method is analyzed. It is shown that the error
from the truncated LG quadrature is independent of the stepsize, and can be
made arbitrarily small by choosing suitable parameters that are given
explicitly. Numerical examples are presented to verify the effectiveness of the
current fast method.
