How to calculate the error of the series that originates the Euler's constant

I have been asked to calculate how many terms of the series that defines the Euler's constsant $\gamma$ to add at least, to calculate the value of $\gamma$ with error less or equal to $3 \times 10^{-3}$.

I have tried several approaches knowing that the error between the value of a series and a partial addition is defined by

$$\left| \sum_{k=1}^{\infty}{f(k)} - \sum_{k=1}^{n}{f(k)} \right| = \left| S-S_n \right| \leq \int_{n}^{\infty}{f(t)dt}$$

and just substituting the definition of $\gamma$ as

$$\gamma = \lim_{n \rightarrow \infty}{\left( \sum_{k=1}^{n}{\dfrac{1}{k} - \log{n}} \right)}$$

but I don't know how to calculate the value of $n$ from there... Any help or suggestion of solving this in other way will be gratefully welcome.

There is the the asymptotic estimation $$H_n = \log n + \gamma + \frac1{2n} - \frac1{2n^2} + \frac1{120n^4} + O\left(\frac1{n^6}\right)$$ (see Introduction to Analytic Number Theory by Apostol)
$$H_n = \log n + \gamma + \frac1{2n} - \sum_{k=1}^\infty\frac{B_{2k}}{2kn^{2k}}$$ with $B_k$ the Bernoulli numbers.