  We present an efficient algorithm for the evaluation of the Caputo fractional
derivative $_0^C\!D_t^\alpha f(t)$ of order $\alpha\in (0,1)$, which can be
expressed as a convolution of $f'(t)$ with the kernel $t^{-\alpha}$. The
algorithm is based on an efficient sum-of-exponentials approximation for the
kernel $t^{-1-\alpha}$ on the interval $[\Delta t, T]$ with a uniform absolute
error $\varepsilon$, where the number of exponentials $N_{\text{exp}}$ needed
is of the order
  $O\left(\log\frac{1}{\varepsilon}\left(
  \log\log\frac{1}{\varepsilon}+\log\frac{T}{\Delta t}\right)
  +\log\frac{1}{\Delta t}\left(
  \log\log\frac{1}{\varepsilon}+\log\frac{1}{\Delta t}\right)
  \right)$.
  As compared with the direct method, the resulting algorithm reduces the
storage requirement from $O(N_T)$ to $O(N_{\text{exp}})$ and the overall
computational cost from $O(N_T^2)$ to $O(N_TN_{\text{exp}})$ with $N_T$ the
total number of time steps. Furthermore, when the fast evaluation scheme of the
Caputo derivative is applied to solve the fractional diffusion equations, the
resulting algorithm requires only $O(N_SN_{\text{exp}})$ storage and
$O(N_SN_TN_{\text{exp}})$ work with $N_S$ the total number of points in space;
whereas the direct methods require $O(N_SN_T$) storage and $O(N_SN_T^2)$ work.
The complexity of both algorithms is nearly optimal since $N_{\text{exp}}$ is
of the order $O(\log N_T)$ for $T\gg 1$ or $O(\log^2N_T)$ for $T\approx 1$ for
fixed accuracy $\varepsilon$. We also present a detailed stability and error
analysis of the new scheme for solving linear fractional diffusion equations.
The performance of the new algorithm is illustrated via several numerical
examples. Finally, the algorithm can be parallelized in a straightforward
manner.
