Prove closed expression for $\int_0^1 \log(x) \log(1+x) \log(2+x)\,dx$ In Closed form of $\int_{0}^{1} \frac{\log(1+x)\log(2+x) \log(3+x)}{1+x}\,dx$ I have proposed an integral which I could not solve, and though there were some upvotes on the question no solution was provided. Hence I looked for simplifications which still are not trivial. 
Here's an example where I found a closed expression with the help of Mathematica which can be verified numerically but I'm lacking a proof. 
Hence my question is
Prove that
$$\int_0^1 \log(x)\log(x+1)\log(x+2)\,dx \\
= -6+3 \log ^3(2)-\frac{\log ^3(3)}{3}+\frac{\log ^2(2)}{2}-3 \log (3) \log (2)+6 \log (3)\\+\zeta(2) (1-2 \log (2))-\frac{13 \zeta (3)}{8}\\-\operatorname{Li}_2\left(-\frac{1}{2}\right)-6 \operatorname{Li}_2\left(-\frac{1}{2}\right) \log (2)+4 \operatorname{Li}_2\left(\frac{1}{4}\right) \log (2)\\-2 \operatorname{Li}_2\left(\frac{1}{3}\right) \log (3)+\operatorname{Li}_2\left(-\frac{1}{3}\right) \log (3)\\
-4 \operatorname{Li}_3\left(-\frac{1}{2}\right)-2 \operatorname{Li}_3\left(\frac{1}{3}\right)+\operatorname{Li}_3\left(-\frac{1}{3}\right)+2 \operatorname{Li}_3\left(\frac{1}{4}\right)\\\simeq -0.18403235664237885896 $$
Notice that the expression is composed of $\pi$, $\log(s)$, $\zeta(s)$, and $\operatorname{Li}_{s}(t)$.
Remark 1: Mathematica was able to find the antiderivative but it turned out to contain complex valued summands. These cancelled out numerically but I could not prove mathematically that their contribution vanishes. 
Remark 2: I have not found the present class of integrals (product of logs with successively shifted arguments) in the 60 problems the of the book "(Almost) Impossible Integrals, Sums, and Series" by
Cornel Ioan Valean (https://it.b-ok2.org/book/4996918/0df734) which is famous and frequently quoted in this forum. So this type of problem seems to be new.
 A: We have to calculate the integral
$$i = \int_0^1 \log(x)\log(1+x)\log(2+x)\,dx\tag{1}$$
1. My derivation of the closed expression
First I tried to find the indefinite integral (the antiderivative of the integrand)
$$a(x)=\int \log(x)\log(1+x)\log(2+x)\,dx\tag{2}$$
I was lucky, Mathematica quickly returned an expression the drivative of which gave back the integrand.
It turned out that $a(0)=0$ so that $i = a(1)$. The result is also numerically correct to a good approximation.
The expression $a(1)$ formally still contained an imaginary part.
But this imaginary part turns out to be numerically zero, i.e.
$$a_i = -2 \operatorname{Li}_2\left(\frac{1}{3}\right)+\text{Li}_2\left(-\frac{1}{3}\right)+\frac{\pi ^2}{6}-\frac{1}{2} \log ^2(3)= 0\tag{3}$$
I'm sure that $(3)$ holds exactly but I haven't yet found the dilog relation to prove it.
Notice that this derivation is a valid proof: we have used a heuristic tool to find a solution which could be verified.
2. attempt using parametric derivatives, double series
My first solution attempt starts with generating the logs by differentiating the function
$$f=x^a (x+1)^b (x+2)^c$$
with respect to the parameters $a$, $b$, and $c$, and then letting the parameters go to $0$.
Let us expand $f$ into a double binomial series
$$f_s = 2^c x^a \sum _{n=0}^{\infty } \sum _{m=0}^{\infty } \frac{x^m x^n \binom{b}{m} \binom{c}{n}}{2^n}$$
performing the integral gives for the summand
$$s(n,m)=\frac{2^{c-n} \binom{b}{m} \binom{c}{n}}{a+m+n+1}$$
The drivatives and the respective limits are
$$s_a=\frac{\partial s(n,m)}{\partial a}|_{a\to 0} = -\frac{2^{c-n} \binom{b}{m} \binom{c}{n}}{(m+n+1)^2}$$
$$s_b = \frac{\partial s_a}{\partial b}|_{b\to 0} = -\frac{\binom{0}{m} 2^{c-n} (-\psi ^{(0)}(1-m)-\gamma ) \binom{c}{n}}{(m+n+1)^2}$$
$$s_c = \frac{\partial s_b}{\partial c}|_{c\to 0} =-\frac{2^{-n} \binom{0}{m} \binom{0}{n} H_{-m} \left(H_{-n}-\log (2)\right)}{(m+n+1)^2}$$
We observe that harmonicnumbers have een generated but in a peculiar combination with the binomial coefficient.
We know that $H_{z}$ has simple poles at negatives integer $z$. On the other hand $\binom{0}{k}=0$ at natural $k$. In fact there is cancellation described by the formula
$$\lim_{m\to 0} \, \binom{0}{m} H_{-m}= 0$$
$$\lim_{m\to 1} \, \binom{0}{m} H_{-m}=\frac{(-1)^m}{m}$$
For $n=0$ the summand becomes
$$\lim_{n\to 0} \, -\frac{2^{-n} \binom{0}{m} \binom{0}{n} H_{-m} \left(H_{-n}-\log (2)\right)}{(m+n+1)^2}=\frac{\log (2) \binom{0}{m} H_{-m}}{(m+1)^2}$$
So that the remaining $m$-sum starts at $m=1$ and gives
$$\sum _{m=1}^{\infty } \frac{(-1)^m \log (2)}{m (m+1)^2}=\left(-\frac{\pi ^2}{12}+2-2 \log (2)\right) \log (2) $$
Now the true double sum has $n\ge1$, $m\ge1$ so that $\log (2) \binom{0}{n}=0$ and the sum becomes
$$-\sum _{n=1}^{\infty } \sum _{m=1}^{\infty } \frac{2^{-n} (-1)^{m+n}}{m n (m+n+1)^2}$$
I'm just seeing that I have made a simple thing compilcated. We better expand the two logs with the shift into a power series ...
(to be continued).
A: Incomplete solution
First write $\ln(2+x)=\ln2+\ln(1+x/2)$
$$\Longrightarrow I=\ln2\int_0^1 \ln x\ln(1+x)\ dx+\int_0^1\ln x\ln(1+x)\ln(1+x/2)\ dx$$
$$=\ln2 I_1+I_2$$
Apply integration by parts for $I_1$
$$I_1=(x\ln x-x)\ln(1+x)|_0^1-\int_0^1\frac{x\ln x-x}{1+x}\ dx=\boxed{2-2\ln2-\frac12\zeta(2)}$$
For $I_2$, write $\displaystyle\ln(1+x/2)=-\sum_{n=1}^\infty\frac{(-1)^n}{2^nn}x^n=\frac12\sum_{n=2}^\infty\frac{(-1)^n}{2^n(n-1)}x^{n-1}$
$$\Longrightarrow I_2=\frac12\sum_{n=2}^\infty\frac{(-1)^n}{2^n(n-1)}\int_0^1 x^{n-1}\ln x\ln(1+x)\ dx$$
Now use the identity
$$\int_0^1\frac{x^{2n}}{1+x}dx=\ln2+H_n-H_{2n}$$
By integration by parts we have
$$\int_0^1 x^{2n-1}\ln(1+x)\ dx=\frac{H_{2n}-H_n}{2n}$$
replace $2n$ with $n$ then differentiate with respect to $n$
$$\int_0^1 x^{n-1}\ln x\ln(1+x)\ dx=\frac{H_{n/2}}{n^2}+\frac{H_{n/2}^{(2)}}{2n}-\frac{H_n}{n^2}-\frac{H_n^{(2)}}{n}+\frac{\zeta(2)}{2n}$$
Therefore 
$$I_2=\frac12\sum_{n=2}^\infty\frac{(-1)^nH_{n/2}}{2^n(n-1)n^2}+\frac14\sum_{n=2}^\infty\frac{(-1)^nH_{n/2}^{(2)}}{2^n(n-1)n}-\frac12\sum_{n=2}^\infty\frac{(-1)^nH_{n}}{2^n(n-1)n^2}-\frac12\sum_{n=2}^\infty\frac{(-1)^nH_n^{(2)}}{2^n(n-1)n}\\
+\frac{\zeta(2)}{4}\sum_{n=2}^\infty\frac{(-1)^n}{2^n(n-1)n}$$
For the first and second sum, we can use 
$$\sum_{n=2}^\infty f(n)=\sum_{n=1}^\infty f(2n)+\sum_{n=1}^\infty f(2n+1)$$
