We divide the solution into several steps.

1. Reduction to Euler series.

The key ingredient for the reduction is the following integral.

\begin{equation*} \int_{0}^{\frac{\pi}{2}} \sin^{p} \theta \cos^{q} \theta \log \sin\theta \log \cos\theta \; d\theta. \end{equation*} To find a closed form of it, we consider the beta function identity. \begin{equation*} 2 \int_{0}^{\frac{\pi}{2}} \sin^{2z-1} \theta \cos^{2w-1} \theta \; d\theta = \beta(z, w) = \frac{\Gamma(z)\Gamma(w)}{\Gamma(z+w)} \end{equation*} Differentiating with respect to $z$, we obtain \begin{equation*} 4 \int_{0}^{\frac{\pi}{2}} \sin^{2z-1} \theta \cos^{2w-1} \theta \log \sin\theta \; d\theta = \beta(z, w) \left\{ \psi_0 (z) - \psi_0 (z+w) \right\}, \end{equation*} where $\psi_k$ denotes the $k$-th polygamma function. Differentiating with respect to $w$, we finally obtain \begin{align} & 8 \int_{0}^{\frac{\pi}{2}} \sin^{2z-1} \theta \cos^{2w-1} \theta \log \sin\theta \log\cos\theta \; d\theta \nonumber \\ & \qquad = \beta(z, w) \left[ \left\{ \psi_0 (z) - \psi_0 (z+w) \right\}\left\{ \psi_0 (w) - \psi_0 (z+w) \right\} - \psi_1 (z+w) \right], \label{eqn:key_int} \end{align} the very closed form for \eqref{eqn:key_int}.

With this formula we want to reduce \eqref{eqn:wts} to a formula involving a certain type of infinite summations, namely the alternating Euler series. Our key formula \eqref{eqn:key_int} yields \begin{eqnarray} & & \int_{0}^{1} \log (1-x) \log x \log (1+x) \; dx \\ & = & \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n} \int_{0}^{1} x^n \log (1-x) \log x \; dx \\ & = & \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n} \int_{0}^{\frac{\pi}{2}} 8 \sin^{2n+1} \theta \cos \theta \log \sin\theta \log \cos\theta \; dx \qquad (x = \sin^2 \theta) \\ & = & \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n} \beta(n+1, 1) \left[ \{ \psi_0 (n+1) - \psi_0 (n+2) \} \{ \psi_0 (1) - \psi_0 (n+2) \} - \psi_1 (n+1) \right] \\ & = & \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n(n+1)} \left[ \frac{H_{n+1}}{n+1} - \zeta(2) + H_{n+1}^{(2)} \right], \end{eqnarray} where \begin{equation*} H_n = \sum_{k=1}^{n} \frac{1}{k} \quad \text{and} \quad H_n^{(s)} = \sum_{k=1}^{n} \frac{1}{k^s} \end{equation*} are the $n$-th harmonic number and the $n$-th generalized harmonic number of order $s$, respectively, and we have exploited the following identities involving polygamma functions. \begin{equation*} \psi_0 (n+1) = -\gamma + H_n , \qquad \psi_1 (n+1) = \zeta(2) - H_n^{(2)}. \end{equation*}

Now consider the following family of alternating series \begin{equation*} A_{m}^{(s)} = \sum_{n=1}^{\infty} \frac{(-1)^{n-1} H_n^{(s)} }{n^m}, \end{equation*} namely the alternating Euler series. By some algebra, we obtain \begin{align*} \frac{(-1)^{n-1}}{n(n+1)} \left[ \frac{H_{n+1}}{n+1} - \zeta(2) + H_{n+1}^{(2)} \right] & = (-1)^{n-1} \Bigg[ \frac{H_n}{n} - \frac{H_{n+1}}{n+1} - \frac{H_{n+1}}{(n+1)^2} + \frac{H_n^{(2)}}{n} - \frac{H_{n+1}^{(2)}}{n+1} \\ & \qquad \qquad + \left( \frac{1}{n} - \frac{1}{n+1} \right) (2 - \zeta(2)) - \frac{1}{(n+1)^2} \Bigg], \end{align*} from which the following reduction formula follows. \begin{align} &\int_{0}^{1} \log (1-x) \log x \log (1+x) \; dx \nonumber \\ &\qquad = 2A_1^{(1)} + 2A_1^{(2)} + A_2^{(1)} + \zeta(2) \left( \frac{3}{2} - 2\log 2 \right) + 4 \log 2 - 6 \label{eqn:reduce} \end{align} Therefore it suffices to evaluate the following 3 series \begin{equation*} A_1^{(1)}, \quad A_2^{(1)}, \quad \text{and} \quad A_1^{(2)} \end{equation*} in the identity \eqref{eqn:reduce}.

2. Evaluation of $A_m^{(s)}$.

For the evaluation, we make the following observation.

Consider $a_n = \frac{1}{n}$ with $a_0 = 0$. Then \begin{equation*} -\frac{\log(1-x)}{1-x} = \sum_{n=1}^{\infty} H_n x^n. \end{equation*} Integrating this identity, \begin{align*} \sum_{n=1}^{\infty} \frac{H_n}{n} x^n & = - \int_{0}^{x} \frac{\log(1-t)}{t(1-t)} \; dt = - \int_{0}^{x} \left( \frac{\log(1-t)}{t} + \frac{\log(1-t)}{1-t} \right) \; dt \\ & = \left[ \mathrm{Li}_2 (t) + \frac{1}{2} \log^2 (1-t) \right]_{0}^{x} = \mathrm{Li}_2 (x) + \frac{1}{2} \log^2 (1-x). \end{align*} Finally, plugging $x = -1$ yields \begin{equation} A_{1}^{(1)} = \frac{1}{2}\zeta(2) - \frac{1}{2}\log^2 2. \label{eqn:iden_1} \end{equation}

Let's review on the idea in the calculation of \eqref{eqn:iden_1}. We first identified the series as a value of a power series evaluated at some point, and then derived the general formula for the power series. This small triumph seems to suggest that our strategy would also work for other two summations. However, \begin{align} \sum_{n=1}^{\infty} \frac{H_n}{n^2} x^n & = \int_{0}^{x} \frac{1}{t} \sum_{n=1}^{\infty} \frac{H_n}{n} t^n \; dt = \int_{0}^{x} \frac{1}{t} \left( \mathrm{Li}_2 (t) + \frac{1}{2} \log^2 (1-t) \right) \; dt \nonumber \\ & = \mathrm{Li}_3 (x) + \frac{1}{2} \int_{0}^{x} \frac{\log^2 (1-t)}{t} \; dt \label{eqn:iden_2} \end{align} and also \begin{align} \sum_{n=1}^{\infty} \frac{H_n^{(2)}}{n} x^n & = \int_{0}^{x} \frac{1}{t} \sum_{n=1}^{\infty} H_n^{(2)} t^n \; dt = \int_{0}^{x} \frac{1}{t} \left( \frac{\mathrm{Li}_2 (t)}{1-t} \right) \; dt \nonumber \\ & = \int_{0}^{x} \left( \frac{\mathrm{Li}_2 (t)}{t} + \frac{\mathrm{Li}_2 (t)}{1-t} \right) \; dt = \mathrm{Li}_3 (x) + \int_{0}^{x} \frac{\mathrm{Li}_2 (t)}{1-t} \; dt \nonumber \\ & = \mathrm{Li}_3 (x) + \left[ - \mathrm{Li}_2 (t) \log (1 - t) \right]_{0}^{x} - \int_{0}^{x} \frac{\log^2 (1 - t)}{t} \; dt \nonumber \\ & = \mathrm{Li}_3 (x) - \mathrm{Li}_2 (x) \log (1 - x) - \int_{0}^{x} \frac{\log^2 (1 - t)}{t} \; dt \label{eqn:iden_3} \end{align} shows that this idea is thwarted by the same factor \begin{equation}\label{eqn:int_in_q} I(x) = \frac{1}{2} \int_{0}^{x} \frac{\log^2 (1-t)}{t} \; dt. \end{equation} Therefore it is essential for our calculation to calculate the value of \eqref{eqn:int_in_q} at $x = -1$.

In the calculation \eqref{eqn:iden_3}, we actually derived the following identity \begin{equation*} \int_{0}^{x} \frac{\mathrm{Li}_2 (t)}{1-t} \; dt \ = \ - \mathrm{Li}_2 (x) \log (1 - x) - 2I(x). \end{equation*} Now we are going to evaluate the left-hand side of this identity in another way. Using Euler's reflection formula for the dilogarithm, \begin{align*} \int_{0}^{x} \frac{\mathrm{Li}_2 (t)}{1-t} \; dt & = \int_{0}^{x} \frac{1}{1-t} \left[ \zeta(2) - \log t \log (1-t) - \mathrm{Li}_2 (1-t) \right] \; dt \\ & = \int_{0}^{x} \frac{1}{1-t} \left[ \zeta(2) - \mathrm{Li}_2 (1-t) \right] \; dt - \int_{0}^{x} \frac{\log t \log (1-t)}{1-t} \; dt \\ & = \int_{1-x}^{1} \frac{\zeta(2) - \mathrm{Li}_2 (u)}{u} \; du - \int_{0}^{x} \frac{\log t \log (1-t)}{1-t} \; dt \qquad (u = 1-t) \\ & = - \zeta(2) \log (1-x) - \zeta(3) + \mathrm{Li}_3 (1-x) - \int_{0}^{x} \frac{\log t \log (1-t)}{1-t} \; dt \\ & = - \zeta(2) \log (1-x) - \zeta(3) + \mathrm{Li}_3 (1-x) + \left[ \frac{1}{2} \log^2 (1-t) \log t \right]_{0}^{x} - I(x) \\ & = - \zeta(2) \log (1-x) - \zeta(3) + \mathrm{Li}_3 (1-x) + \frac{1}{2} \log^2 (1-x) \log x - I(x). \end{align*} From this calculation, we obtain a closed for for \eqref{eqn:int_in_q} as follows: \begin{equation}\label{eqn:int_in_q2} I(x) \ = \ \zeta(3) - \mathrm{Li}_3 (1-x) + [\zeta(2) - \mathrm{Li}_2 (x)] \log (1 - x) - \frac{1}{2} \log^2 (1-x) \log x. \end{equation} This identity, however, does not immediately yield the value for \eqref{eqn:iden_2} and \eqref{eqn:iden_3}. We have to evaluate the value of $I(x)$ at $x = -1$, but there are some terms in \eqref{eqn:int_in_q2} which are undefined at $x = -1$. This means that, in order to utilize this formula, we have to analyze the behavior of each terms of \eqref{eqn:int_in_q2} more carefully.

3. Principal Value of the Polylogarithm.

At the core of the problem lies the branch cut, namely $(-\infty, 0)$, of $\mathrm{Li}_3 (1-x)$ and $\log x$. To circumvent this problem, we consider $I(x)$ as a function on a region in the complex plane. Since the branch cut of the integrand is $[1, \infty)$ and $I(1)$ is defined, \eqref{eqn:int_in_q} defines a holomorphic function on $\mathbb{C}-(1,\infty)$. By analytic continuation, the identity \eqref{eqn:int_in_q2} remains valid whenever both sides are defined. Therefore, instead of plugging $x = -1$ to \eqref{eqn:int_in_q2} directly, we can evaluate $I(-1)$ by taking limit to \eqref{eqn:int_in_q2} along a certain path.

Here, we choose the path as $x = -1 - ih$ as $h \to 0^+$. From Euler's reflection formula again, \begin{align*} \mathrm{Li}_2 (2+0^+ i) & = \mathrm{Li}_2 (1 - (-1 - 0^+ i)) \\ & = \zeta(2) - \log (-1 - 0^+ i) \log (1 - (-1 - 0^+ i)) - \mathrm{Li}_2 (-1 - 0^+ i) \\ & = \zeta(2) - (-i \pi) \log 2 - \text{Li}_2 (-1) \\ & = \frac{3}{2} \zeta(2) + i \pi \log 2 . \end{align*} Also, plugging $x = -1-0^+ i$ and $x = \frac{1}{2}$ respectively to the following trilogarithm identity \begin{equation*} \mathrm{Li}_3 (x) + \mathrm{Li}_3 (1-x) + \mathrm{Li}_3 \big( 1 - \tfrac{1}{x} \big) = \zeta(3) + \zeta(2) \log x + \tfrac{1}{6}\log^3 x - \tfrac{1}{2}\log^2 x \log(1-x), \end{equation*} we obtain \begin{equation*} \text{Li}_3 (2+0^+ i) \ = \ \frac{7}{8}\zeta(3) + \frac{3}{2}\zeta(2) + \frac{i\pi}{2} \log^2 2. \end{equation*} Therefore, plugging $x = -1-0^+ i$ to \eqref{eqn:int_in_q2} and simplifying, we obtain \begin{equation*} I(-1) = \frac{1}{8} \zeta(3) \end{equation*} Therefore, from \eqref{eqn:iden_2} and \eqref{eqn:iden_3}, we obtain \begin{equation*} A_2^{(1)} \ = \ \frac{5}{8} \zeta(3) , \qquad A_1^{(2)} \ = \ \zeta(3) - \frac{1}{2}\zeta(2)\log 2. \end{equation*}

4. Conclusion.

Summarizing all our efforts, we obtain \begin{align*} \int_{0}^{1} \log (1-x) \log x \log (1+x) \; dx & = \frac{21}{8} \zeta(3) + \frac{1}{2} \zeta(2) (5 - 6 \log 2) + 4 \log 2 - \log^2 2 - 6. \end{align*}

References

1. Integrals in "AoPS Forum"
2. An Interesting Integral in "AoPS Forum"
3. Philippe Flajolet and Bruno Salvy, "Euler Sums and Contour Integral Representations".