So one of the approaches to proving the equality in the question is via the following three steps: First differentiate both sides of the equation to see that they agree up to a constant. This reduces to showing the case of $x = 1$, for which $\log x = 0$. Next we apply integration by parts to get \begin{align*} \int_1^{\infty} \exp(-y)/y dy - \int_0^1 \frac{1-\exp(-y)}{y}dy = \int_0^{\infty} \exp(-y) \log y dy \end{align*} Finally observe that $\Gamma'(1)$ equals the RHS, by differentiating under the integral sign, valid because things are decaying fast enough at infinity. So it remains to show $\Gamma'(1) =\gamma$. I saw a soft argument (i.e., without using infinite product) in the link scipp.ucsc.edu/~haber/ph116A/psifun_10.pdf This is re-exposed below: first we establish that for $\Psi(x) = \log \Gamma(x)$, \begin{align*} \Psi'(x+1) = \Psi'(x) + 1/x \end{align*} This is easy enough since we have we have the functional equation $\Gamma(x+1) = x\Gamma(x)$. Next using stirling approximation we get \begin{align*} \Psi(x+1) = (x+1/2)\log x -x + 1/2 \log 2 \pi + O(1/x) \end{align*} and then they differentiate this and claim that $O(1/x)' = O(1/x^2)$, which is clearly false (take $f(x) = 1/x cos(e^x)$). But I found in Wikipedia another formula that gives the precise error term in terms of an integral of the monotone function $arctan(1/x)$. So this is enough to establish $O(1/x^2)$ for the error term in the derivative of $\Psi$. So we get the asymptotics $\lim_{x \to \infty} \Psi'(x+1) = \log(x)$, from which we get $\Psi'(1) = \gamma$. Now notice $\Psi'(x) = \Gamma'(x)/ \Gamma(x)$, and $\Gamma(1) = 1$, so $\Gamma'(1) = \gamma$ also.