  We provide efficient methods to evaluate the Riemann zeta, the Lerch zeta and
the Dirichlet $L$-functions. The method uses the Riemann-Siegel (RS) type
formulas and a modified double exponential (MDE) quadrature method near the
saddle point of appropriate integrands. We provide simplified derivations of
the RS formulas, containing finite series sums and residual integrals, for the
Lerch and the Dirichlet $L$-functions. The MDE method allows us to remove the
contribution of the singularities near the contour of integration for the
residual integrals. The method allows us to evaluate the residual integrals to
any prescribed accuracy. The numerical cost of evaluating them is minimal
compared to the main series sums. Thus even highly oscillatory integrands for
these functions can be evaluated with the same complexity as the RS formula. In
particular, the numerical complexity to evaluate the $\zeta(s)$ and the Lerch
zeta function with $s=\sigma+i t$ scales as $\sqrt{t}$, and for Dirichlet
$L$-function it scales as $\max(q, \sqrt{q t})$. The method allows the
automatic setting of integral cutoffs and finite discretization size to achieve
prescribed accuracy. Furthermore, it ensures that every halving of the
discretization interval almost leads to doubling the accuracy of the results.
Thus, it allows for a further check on the accuracy of the results obtained.
