We introduce an algorithm to compute the functions belonging to a suitable set F defined as follows: f∈F means that f(s,x), s∈A⊂R being fixed and x>0, has a power series expansion centred at x0=1 with convergence radius greater or equal than 1; moreover, it satisfies a functional equation of step 1 and the Euler-Maclaurin summation formula can be applied to f. Denoting the Euler gamma-function as Γ, we will show that, for x>0, logΓ(x), the digamma function ψ(x), the polygamma functions ψ(w)(x), w∈N, w≥1, and, for s>1 being fixed, the Hurwitz ζ(s,x)-function and its first partial derivative ∂ζ∂s(s,x) are in F. In all these cases the coefficients of the involved power series will depend on the values of ζ(u), u>1, where ζ is the Riemann zeta-function. As a by-product, we will also show how to compute the Dirichlet L-functions L(s,χ) and L′(s,χ), s>1, χ being a primitive Dirichlet character, by inserting the reflection formulae of ζ(s,x) and ∂ζ∂s(s,x) into the first step of the Fast Fourier Transform algorithm. Moreover, we will obtain some new formulae and algorithms for the Dirichlet β-function and for the Catalan constant G. Finally, we will study the case of the Bateman G-function and of the alternating Hurwitz zeta-function, also known as the η-function; we will show that, even if they are not in F, our approach can be adapted to handle them too. In the last section we will also describe some tests that show a performance gain with respect to a standard multiprecision implementation of ζ(s,x) and ∂ζ∂s(s,x), s>1, x>0.