  We express the Riemann zeta function $\zeta\left(s\right)$ of argument
$s=\sigma+i\tau$ with imaginary part $\tau$ in terms of three absolutely
convergent series. The resulting simple algorithm allows to compute, to
arbitrary precision, $\zeta\left(s\right)$ and its derivatives using at most
$C\left(\epsilon\right)\left|\tau\right|^{\frac{1}{2}+\epsilon}$ summands for
any $\epsilon>0$, with explicit error bounds. It can be regarded as a
quantitative version of the approximate functional equation. The numerical
implementation is straightforward. The approach works for any type of zeta
function with a similar functional equation such as Dirichlet $L$-functions, or
the Davenport-Heilbronn type zeta functions.
