Generating function of square harmonic numbers During my study of generating functions, I was able to calculate the generating function of the sequence of harmonic numbers $H_n$:
$$\sum_{n=1}^\infty H_nx^n=\frac{\ln(1-x)}{x-1}$$
However, I also tried to find generating functions for $H_n^2$ and $H_n^3$, with which I was unsuccessful (the rearrangement method I used for the generating function of $H_n$ didn't reduce as nicely for $H_n^2$ and $H_n^3$). 
Any hints about how to find
$$\sum_{n=1}^\infty H_n^2x^n=\space ?$$
and
$$\sum_{n=1}^\infty H_n^3x^n=\space ?$$
Please don't write a full answer and spoil it for me - I just want a hint.
 A: We have that
$$ f(x)=\sum_{n\geq 0}a_n x^n\quad\Longleftrightarrow\quad \frac{f(x)}{1-x}=\sum_{n\geq 0}A_n x^n $$
where $A_n = a_0+a_1+\ldots+a_n$. In order to find the OGF of $H_{n}^2$ it is enough to find the OGF of
$$ H_{n+1}^2-H_{n}^2 = \left(H_{n+1}-H_n\right)\left(H_{n+1}+H_n\right)=\frac{2H_n}{n+1}+\frac{1}{(n+1)^2}, $$
and while the OGF of $\frac{1}{(n+1)^2}$ is clearly related with $\text{Li}_2(x)=\sum_{n\geq 1}\frac{x^n}{n^2}$, the OGF of $\frac{H_n}{n+1}$ can be deduced by applying termwise integration to the OGF of $H_n$. It follows that
$$ \sum_{n\geq 1}H_n^2 x^n = \frac{\log^2(1-x)+\text{Li}_2(x)}{1-x}\tag{A}$$
and the (more involved) OGF of $H_n^3$ can be computed through the same trick.
It is useful to consider that the Taylor series of $\log(1-x)^k$ is related to Stirling numbers of the first kind.
A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

Hereafter, $\ds{\bracks{\cdots}}$ is an Iverson Bracket. Namely, $\ds{\bracks{P} = 1}$ whenever $\ds{P\ \mbox{is}\ \color{red}{\texttt{true}}}$ and $\ds{0\  \color{red}{\mbox{otherwise}}}$.

\begin{align}
\left.\sum_{n = 1}^{\infty}H_{n}^{2}x^{n}\,\right\vert_{\ \verts{x}\ <\ 1} & =
\sum_{n = 1}^{\infty}\overbrace{\braces{\sum_{i = 1}^{\infty}{\bracks{i \leq n} \over i}}}^{\ds{H_{n}}}\
\overbrace{\braces{\sum_{j = 1}^{\infty}{\bracks{j \leq n} \over j}}}^{\ds{H_{n}}}\ x^{n}
\\[5mm] & =
\sum_{i = 1}^{\infty}{1 \over i}\sum_{j = 1}^{\infty}{1 \over j}
\sum_{n = 1}^{\infty}\bracks{n \geq i}\bracks{n \geq j}x^{n}
\\[5mm] & =
\sum_{i = 1}^{\infty}{1 \over i}\sum_{j = 1}^{\infty}{1 \over j}\braces{%
\bracks{i \leq j}\sum_{n = j}^{\infty}x^{n} +
\bracks{i > j}\sum_{n = i}^{\infty}x^{n}}
\\[5mm] & =
\sum_{i = 1}^{\infty}{1 \over i}\sum_{j = 1}^{\infty}{1 \over j}\braces{%
\bracks{j \geq i}\,{x^{j} \over 1 - x} +
\bracks{j < i}\,{x^{i} \over 1 - x}}
\\[5mm] & =
{1 \over 1 - x}\sum_{i = 1}^{\infty}{1 \over i}\pars{%
\sum_{j = i}^{\infty}{x^{j} \over j} + \sum_{j = 1}^{i - 1}{x^{i} \over j}} =
{1 \over 1 - x}\pars{\sum_{i = 1}^{\infty}{1 \over i}
\sum_{j = i}^{\infty}{x^{j} \over j} +
\sum_{i = 1}^{\infty}{x^{i} \over i}\overbrace{\sum_{j = 1}^{i - 1}{1 \over j}}
^{\ds{H_{i - 1}}}}
\\[5mm] & =
{1 \over 1 - x}\pars{\sum_{j = 1}^{\infty}{x^{j} \over j}\
\overbrace{\sum_{i = 1}^{j}{1 \over i}}^{\ds{H_{j}}}\ +
\sum_{i = 1}^{\infty}{x^{i + 1} \over i + 1}\,H_{i}} =
{1 \over 1 - x}\sum_{i = 1}^{\infty}
\pars{{x^{i} \over i} + {x^{i + 1} \over i + 1}}H_{i}
\\[5mm] & =
{1 \over 1 - x}\sum_{i = 1}^{\infty}H_{i}
\pars{x^{i}\int_{0}^{1}t^{i - 1}\,\dd t + x^{i + 1}\int_{0}^{1}t^{i}\,\dd t}
\\[5mm] & =
{1 \over 1 - x}
\pars{\int_{0}^{1}{1 \over t}\sum_{i = 1}^{\infty}H_{i}\pars{xt}^{i}\,\dd t + x\int_{0}^{1}\sum_{i = 1}^{\infty}H_{i}\pars{xt}^{i}\,\dd t}
\\[5mm] & =
{1 \over 1 - x}\int_{0}^{1}{1 + xt \over t}
\bracks{-\,{\ln\pars{1 - xt} \over 1 - xt}}\,\dd t =
-\,{1 \over 1 - x}\int_{0}^{x}{1 + t \over t\pars{1 - t}}\,
\ln\pars{1 - t}\,\dd t
\end{align}

because $\ds{\sum_{i = 1}^{\infty}H_{i}z^{i} =
-\,{\ln\pars{1 - z} \over 1 - z}:\ \pars{~H_{i}\ Generating\ Function~}}$.

Then,
\begin{align}
\left.\sum_{n = 1}^{\infty}H_{n}^{2}x^{n}\,\right\vert_{\ \verts{x}\ <\ 1} & =
-\,{1 \over 1 - x}\bracks{%
2\
\underbrace{\int_{0}^{x}{\ln\pars{1 - t} \over 1 - t}\,\dd t}
_{\ds{-\,{1 \over 2}\,\ln^{2}\pars{1 - x}}}\ +\
\underbrace{\int_{0}^{x}\overbrace{{\ln\pars{1 - t} \over t}}
^{\ds{-\,\mrm{Li}_{2}'\pars{x}}}\ \,\dd t}
_{\ds{-\,\mrm{Li}_{2}\pars{x}}}}
\\[5mm] & =
\bbx{\ln^{2}\pars{1 - x} + \mrm{Li}_{2}\pars{x} \over 1 - x}
\end{align}
