Proof of an exact formula for $H_n$ The $n$th harmonic number $H_n$ is defined as
$$H_n=\sum_{n\geq k\geq 1}\frac{1}{k}$$
A good approximation for this is
$$H_n\approx \gamma+\log n +\frac{1}{2n}$$
Where $\gamma$ is the Euler-Mascheroni constant.
In the book Table of integrals, infinite series, and products by I.S Gradshteyn and A.M Ryzhik, it is given that
$$H_n=\gamma+\log n+\frac{1}{2n}+\sum_{k\geq 2}\frac{A_n}{n(n+1)...(n+k-1)}$$
where
$$A_k=\frac{1}{k}\int_{0}^{1}x(1-x)(2-x)...(k-1-x)dx$$
This is an exact formula. I found it pretty remarkable. How can one prove this?
There is another question. The starting values of $A_n$ are $0,\frac{1}{12},\frac{1}{12},\frac{9}{20}$. Is there any pattern in these numbers? Another way to say this it: is there any closed form for the integral?
Update: I found in Wolfram mathworld that
$$H_n=\gamma+\psi_{0}(n+1)$$
Where $\psi_{0}(x)$ is the digamma function. I did a little research and found
$$\psi_{0}(1+z)=\log(z)+\frac{1}{2z}-\sum_{j\geq 1}\frac{B_{2j}}{2jz^{2j}}$$
so
$$H_n=\gamma+\log(n)+\frac{1}{2n}-\sum_{j\geq 1}\frac{B_{2j}}{2jz^{2j}}$$
so the only thing that is left to prove is
$$\sum_{j\geq 1}\frac{B_{2j}}{2jn^{2j}}=-\sum_{k\geq 2}\frac{A_n}{n(n+1)...(n+k-1)}$$
where $B_{2n}$ are the Bernoulli numbers.
How can we prove the equality of those two series? From this the original formula for $H_n$ can be proved. You can give the proof in any way you like but the proof of the equality of these two series or the proof with the fact that $H_n=\gamma+\psi_{0}(n+1)$ would be the best.
Update: From the wikipedia page of Gregory coefficients (link to this article was given by Donald Splutterwit) I found some properties of $A_k$ from the references of that wikipedia article. This is the article about $A_k$. They are written as $P_{n+1}(y)$ in this paper.
 A: Note that by the Beta function formula we have $$\int_0^1 t^{k-1}(1-t)^{n-1}\,dt = \frac{\Gamma(k)\Gamma(n)}{\Gamma(n+k)} = \frac{(k-1)!}{n(n+1) \cdots (n+k-1)}.$$
Therefore, the series \begin{align} \sum\limits_{k=1}^{\infty} \frac{A_k}{n(n+1) \cdots (n+k-1)} &= \sum\limits_{k=1}^{\infty} \int_0^1 \frac{A_k}{(k-1)!}t^{k-1}(1-t)^{n-1}\,dt \tag{1} \\&= \sum\limits_{k=1}^{\infty} \int_0^1\int_0^1 (-1)^{k-1}\frac{\Gamma(x+1)}{k!\Gamma(x+1-k)}t^{k-1}(1-t)^{n-1}\,dx\,dt \tag{2} \\&= \int_0^1\int_0^1 \left[\sum\limits_{k=1}^{\infty}(-1)^{k-1}\binom{x}{k}t^{k-1}\right](1-t)^{n-1}\,dx\,dt \tag{3} \\&= \int_0^1\int_0^1 \left[\frac{1 - (1-t)^x}{t}\right](1-t)^{n-1}\,dx\,dt \tag{4} \\&= \int_0^1\int_0^1 \left[\frac{1 - t^x}{1-t}\right]t^{n-1}\,dx\,dt \tag{5} \\&= \int_0^1 \left[\frac{1}{1-t} + \frac{1}{\ln t}\right]t^{n-1}\,dt \tag{6} \end{align} where, in line $(2)$ we used the expression for Gamma function and in line $(3)$ we used the generalized Binomial theorem.
The integral $(6)$ is $\displaystyle \int_0^1 \left[\frac{1}{1-t} + \frac{1}{\ln t}\right]t^{n-1}\,dt = H_n - \ln n - \gamma$, a standard Cauchy-Frullani integral.
(will add further details if required! :) )
