Formula for the harmonic series $H_n = \sum_{k=1}^n 1/k$ due to Gregorio Fontana.

My question was inspired by this stackexchange question. For the last 90 minutes I have been trying to prove this formula due to Gregorio Fontana:

$$H_n = \gamma + \log n + {1 \over 2n} - \sum_{k=2}^\infty { (k-1)! C_k \over n(n+1)\ldots(n+k-1)}, \qquad \textrm{ for } n=1,2,3,\ldots,$$

where $H_n = \sum\limits_{k=1}^n 1/k$ and the coefficients $C_k$ are the Gregory coefficients given by $${ z \over \log(1-z)} = \sum_{n=0}^\infty C_k z^k \qquad \textrm{ for } |z|<1.$$

It's a bit frustrating as it's something I recall proving as a student many years ago. I have a vague recollection that I began with something like:

$$H_n = \int_0^1 {1-(1-x)^n \over x } \textrm{d}x,$$

but my attempts to follow on from there have failed. Can you help?

Let's rewrite the formula, using $C_1=1/2$ as $$S_n=\gamma+\log n-H_{n-1}$$ where $$S_n=\sum_{k=1}^\infty\frac{(k-1)!C_k}{n(n+1)\cdots(n+k-1)}.$$ Use the formula $$\frac{(k-1)!}{n(n+1)\cdots(n+k-1)}=\int_0^1 x^{n-1}(1-x)^{k-1}dx$$ which is a special case of the beta integral. Therefore $$S_n=\int_0^1x^{n-1}\sum_{k=1}^\infty C_k(1-x)^k dx=\int_0^1\left(\frac{1}{1-x}+\frac{1}{\log x}\right)x^{n-1}dx.$$ When $n=1$ we get $$S_1=\int_0^1\left(\frac{1}{1-x}+\frac{1}{\log x}\right)dx=\gamma$$ by a well-known formula which can be derived from the aymptotic $$\int_t^\infty\frac{e^{-x}}{x}dx=-\log t-\gamma+O(t)$$ as $t\to0$, for the exponential integral.

By induction it suffices to consider the difference $$S_n-S_{n+1}=\int_0^1\left(x^{n-1}+\frac{x^{n-1}-x^n}{\log x}\right)dx =\frac1n-\int_0^\infty\frac{e^{-ny}-e^{-(n+1)y}}{y}dy.$$ Using the identity $$\int_0^\infty\frac{e^{-ay}-e^{-by}}{y}dy=\log\frac ba$$ for $0 < a < b$ which follows by integrating $e^{-xy}$ over the region $[a,b]\times[0,\infty)$ we get $$S_n-S_{n+1}=\frac1n-\log\frac{n+1}{n}.$$ The desired formula for $S_n$ now follows by induction.