Knuth's article can be found here. Essentially the idea is to compare the approximation for finite $n$ with the equality obtained for $n$ infinite (at this point I want to take issue with the use of "$=$" in the book, which is obviously nonsense since $\gamma$ is constant and the right-hand side clearly has different values for differing values of $N$).
So, finite $n$:
$$ \sum_{s=1}^{n} \frac{1}{s} = \log{n} + \frac{1}{2} + \frac{1}{2n} + \sum_{k=1}^N \frac{B_{2k}}{2k}\left( 1 - \frac{1}{n^{2k}} \right) - \int_1^n \frac{B_{2k+1}(\{x\})}{x^{2k+2}} \, dx \tag{1} $$
(Knuth has written the remainder in a format that makes it easier to deal with here, using the periodic versions of the Bernoulli polynomials. Since it's essentially an integration by parts and a sum, it really is equality. It's going to disappear anyway, so don't worry too much about the exact form.) Taking $n\to \infty$, (1) becomes
$$ \gamma = \frac{1}{2} + \sum_{k=1}^N \frac{B_{2k}}{2k} - \int_1^{\infty} \frac{B_{2k+1}(\{x\})}{x^{2k+2}} \, dx, $$
and subtracting (1) from this gives
$$ \gamma = \sum_{s=1}^{n} \frac{1}{s} - \log{n} - \frac{1}{2n} + \sum_{k=1}^N \frac{B_{2k}}{2k}\frac{1}{n^{2k}} - \int_n^{\infty} \frac{B_{2k+1}(\{x\})}{x^{2k+2}} \, dx. $$
The remainder integral is much smaller than the other terms if we take $n$ a reasonable size (Knuth uses Stirling's formula on an explicit expression for $B_{2k+1}(\{x\})$ to exhibit this), and thus one ends up with the expression in the book.
