Compute an integral containing a product of powers of logarithms. Let $q \ge 1$ and $p \ge 0$ be integers. Consider a following integral:
\begin{equation}
{\mathcal I}^{(q,p)} := \int\limits_0^1 [\log(1-\eta)]^q [\log(\eta)]^p \frac{1}{\eta} d \eta
\end{equation}
Clearly that integral is proportional to the Nielsen generalized polylogarithm at unity.
Now by using Euler's beta function integral it is easy to see that :
\begin{equation}
{\mathcal I}^{(q,p)} =  
\frac{\partial^p}{\partial \theta_1^p} 
\frac{\partial^q}{\partial \theta_2^q} 
\left.
\left[
\frac{(\theta_1-1)! \theta_2!}{(\theta_1+\theta_2)!}
\right]
\right|_{\theta_1=\theta_2=0}
\end{equation}
We have computed the derivative with respect to $\theta_2$ using the Faa di Bruno's formula then we have set $\theta_2=0$ then we differentiated the result $p$ times using Mathematica and finally set $\theta_1=0$. As a result we discovered the following relations:
\begin{eqnarray}
&&1!\cdot{\mathcal I}^{(q,0)} = - \Psi^{(q)}(1) \\
&&2! \cdot{\mathcal I}^{(q,1)} = - \Psi^{(q+1)}(1) + \sum\limits_{j=1}^{q-1} \binom{q}{j} \Psi^{(j)}(1) \Psi^{(q-j)}(1) \\
&&3! \cdot {\mathcal I}^{(q,2)} = -2 \Psi^{(q+2)}(1) + 3 \cdot 1_{q\ge 2}\cdot\sum\limits_{j=1}^{q-1} \binom{q}{j}\cdot\left[
 \Psi^{(j+1)}(1) \Psi^{(q-j)}(1)+\Psi^{(j)}(1) \Psi^{(q+1-j)}(1) \right]+ \\
&&-2 \cdot 1_{q \ge 3}\cdot\sum\limits_{1 \le j < j_1 \le q-1} \binom{q}{j,j_1-j,q-j_1} \Psi^{(j)}(1) \Psi^{(j_1-j)}(1) \Psi^{(q-j_1)}(1)\\
&&4! \cdot {\mathcal I}^{(q,3)} = -6 \Psi^{(q+3)}(1)+\\
&&12 \cdot\sum\limits_{j=1}^{q-1} \binom{q}{j} \left[ \Psi^{(j)}(1) \Psi^{(q-j+2)}(1) + \frac{3}{2} \Psi^{(j+1)}(1) \Psi^{(q-j+1)}(1)+\Psi^{(j+2)}(1) \Psi^{(q-j+0)}(1)\right]+\\
&&-12 \cdot \sum\limits_{1 \le j < j_1 \le q-1} \binom{q}{j,j_1-j,q-j_1} \left[\Psi^{(j)}(1) \Psi^{(j_1-j)}(1) \Psi^{(q-j_1+1)}(1)+\Psi^{(j)}(1) \Psi^{(j_1-j+1)}(1) \Psi^{(q-j_1)}(1)+\Psi^{(j+1)}(1) \Psi^{(j_1-j)}(1) \Psi^{(q-j_1)}(1)\right]+\\
&&6\cdot \sum\limits_{1\le j < j_1 < j_2 \le q-1} \binom{q}{j,j_1-j,j_2-j_1,q-j_2} \Psi^{(j)}(1) \Psi^{(j_1-j)}(1) \Psi^{(j_2-j_1)}(1) \Psi^{(q-j_2)}(1)
\end{eqnarray}
where $\Psi^{(j)}(1)$ is the polygamma function at unity.
Now the question would be how do we find the result for $p \ge 3$? The multitude of terms that appear in the Faa di Bruno formula is difficult to deal with. Is there a more elegant way of arriving at the result? 
 A: We have
$$ \log^q(1-\eta) = q!\sum_{n\geq q}(-1)^q{\,n\, \brack q}\frac{\eta^n}{n!}\tag{1} $$
hence
$$ \mathcal{I}^{(q,p)}=p!q!\sum_{n\geq q}\frac{(-1)^{p+q}}{n!\,n^{p+1}}{\,n\,\brack q}\tag{2} $$
and the problem boils down to computing some Euler sums, once the Stirling numbers of the first kind are converted into combinations of generalized harmonic numbers. In this context values of $p$ or $q$ greater than $3$ lead to intractable problems by hand: that is a good moment for invoking the help of a CAS, without shame.
A: Since understanding of this problem is essential for calculating a whole family of Euler sums I will write down a detailed solution to the problem. We start from rewriting our quantities as follows:
\begin{eqnarray}
{\mathcal I}^{(q,p)} &=& \left.\frac{\partial^p}{\partial \theta_1^p}  e^{\log(\Gamma(\theta_1)} \frac{\partial^q}{\partial \theta_2^q}e^{\log(\Gamma(\theta_2+1)) - \log(\Gamma(\theta_2+\theta_1+1))} \right|_{\theta_1=\theta_2=0}\\
&=& \left.\sum\limits_{1 \cdot m_1+ \cdots + q \cdot m_q = q} \left(\frac{q!}{\prod\limits_{\xi=1}^q m_\xi! (\xi!)^{m_\xi}}\right) \cdot \frac{\partial^p}{\partial \theta_1^p} \left( \frac{\prod\limits_{j=1}^q [\Psi^{(j-1)}(1) - \Psi^{(j-1)}(1+\theta_1)]^{m_j}}{\theta_1} \right) \right|_{\theta_1=0}
\end{eqnarray}
In here we have computed the higher order derivative with respect to $\theta_2$ by using the Faa di Bruno formula. The sum in the right hand side in the equation above runs over all non-negative integers $\left\{ m_j \right\}_{j=1}^q$ such that $\sum\limits_{j=1}^q j m_j = q$. 
Now in the fraction in parentheses under the derivative operator the numerator behaves as $\theta_1^{m_1+m_2+\cdots+m_q}$ as $\theta_1 \rightarrow 0$. This allows us to conclude that only those terms in the sum contribute where $0 \le m_1+m_2+\cdots+m_q \le p+1$.
Bearing this in mind it is not hard to realize that that gigantic sum can be re-written as follows:
\begin{eqnarray}
{\mathcal I}^{(q,p)} = \left.\sum\limits_{l=0}^{p+1} \frac{1}{(l+1)!} \sum\limits_{\begin{array}{r}j_0+j_1+\cdots+j_l=q\\j_0\ge1,\cdots,j_l\ge 1\end{array}} \binom{q}{j_0,\cdots,j_l} \cdot \frac{\partial^p}{\partial \theta^p} \left(\frac{\prod\limits_{\xi=0}^l [ \Psi^{(j_\xi-1)}(1) - \Psi^{(j_\xi-1)}(1+\theta)]}{\theta}\right) \right|_{\theta=0}
\end{eqnarray}
At this stage the only thing we need to do is to evaluate the derivative at zero. It is clear that it is not a hard thing to do. The easiest way to do it is to expand the numerator in parentheses into a Taylor series in $\theta$ and to extract the coefficient at $\theta^{p+1}$. That coefficient multiplied by $p!$ (which comes from differentiating the power function at zero) gives the result. For extreme (i.e. both the smallest and biggest values of $l$ the result has a neat closed form yet for intermediate values of $l$ the number of possible terms that contribute to a given power of $\theta$ increases very rapidly and as such it is difficult to write down the expression in a neat closed form. Bearing all this in mind we actually write down the derivative in question below:
\begin{eqnarray}
&&\left.\frac{\partial^p}{\partial \theta^p} \left( \right) \right|_{\theta=0}=\left\{
\begin{array}{rr}
(-1)^1\frac{1}{(p+1)} \Psi^{(q+p)}(1) & \mbox{if $l=0$}\\
(-1)^2\sum\limits_{l_0+l_1=p+1}\frac{p!}{l_0! l_1!} \prod\limits_{\xi=0}^1 \Psi^{(j_\xi+l_\xi-1)}(1)  & \mbox{if $l=1$}\\
(-1)^3\sum\limits_{l_0+l_1+l_2=p+1}\frac{p!}{l_0! l_1! l_2!} \prod\limits_{\xi=0}^2 \Psi^{(j_\xi+l_\xi-1)}(1)  & \mbox{if $l=2$}\\
\vdots \\
(-1)^{p+1} p! \prod\limits_{\xi=0}^p \Psi^{(j_\xi)}(1) & \mbox{if $l=p$} \\
0 & \mbox{if $l=p+1$}
\end{array}
\right.
\end{eqnarray}
In here the $l$-indices are strictly positive.
This concludes the proof.
Corollary: Consider a slightly more general integral:
\begin{equation}
{\mathcal I}^{(q,p)}_r := \int\limits_0^1 \log(1-\eta)^q \cdot \log(\eta)^p\frac{1}{\eta^{1-r}} d\eta
\end{equation}
where $r\ge 0$. Then by redoing the calculations above we conclude that:
\begin{eqnarray}
{\mathcal I}^{(q,p)}_r= \left.\sum\limits_{1\cdot m_1+\cdots q \cdot m_q = q} \left(\frac{q!}{\prod\limits_{\xi=1}^q m_\xi! (\xi!)^{m_\xi}}\right) \cdot 
\frac{\partial^p}{\partial \theta_1^p} \left( \frac{\prod\limits_{j=1}^q (\Psi^{(j-1)}(1) - \Psi^{(j-1)}(1+\theta_1+r))}{\theta_1+r}\right)\right|_{\theta_1=0}
\end{eqnarray}
The expression on the right hand side is not singular anymore (as it was in the previous case ) and it can be further evaluated by using the chain rule.
