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?
 A: 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.
