Evaluate $\int^1_0 \log^2(1-x) \log^2(x) \, dx$ I have no idea where to even start. WolframAlpha cant compute it either.
$$\int^1_0 \log^2(1-x) \log^2(x) \, dx$$
I think it can be done with series, but I am not sure, can someone help a little so I can get a start??
Thanks!
 A: Using the series of $\ln^2(1-x)$,
\begin{align}
\int^1_0\ln^2{x}\ln^2(1-x) \ {\rm d}x
&=\sum^\infty_{n=1}\frac{2H_n}{n+1}\int^1_0x^{n+1}\ln^2{x}\ {\rm d}x\\
&=\sum^\infty_{n=1}\frac{4H_n}{(n+1)(n+2)^3}
\end{align}
Then integrate $f(z)=\dfrac{(\gamma+\psi_0(-z))^2}{(z+1)(z+2)^3}$ along an infinitely large square. The integral vanishes which implies the sum of its residues is zero. At the positive integers,
\begin{align}
\sum^\infty_{n=1}{\rm Res}(f,n)
&=\sum^\infty_{n=1}\operatorname*{Res}_{z=n}\left[\frac{1}{(z+1)(z+2)^3(z-n)^2}+\frac{2H_n}{(z+1)(z+2)^3(z-n)}\right]\\
&=\sum^\infty_{n=1}\frac{2H_n}{(n+1)(n+2)^3}-\sum^\infty_{n=1}\frac{4n+5}{(n+1)^2(n+2)^4}\\
&=\sum^\infty_{n=1}\frac{2H_n}{(n+1)(n+2)^3}+\frac{\pi^4}{30}+2\zeta(3)-\frac{91}{16}
\end{align}
At $z=0$,
\begin{align}
{\rm Res}(f,0)
&=\operatorname*{Res}_{z=0}\frac{1}{z^2(z+1)(z+2)^3}\\
&=-\frac{5}{16}
\end{align}
At $z=-1$,
\begin{align}
{\rm Res}(f,-1)
&=0
\end{align}
At $z=-2$,
\begin{align}
{\rm Res}(f,-2)
&=\frac{1}{2}\lim_{z\to-2}\frac{{\rm d}^2}{{\rm d}z^2}\frac{(\gamma+\psi_0(-z))^2}{z+1}\\
&=-\frac{\pi^4}{36}+2\zeta(3)+\frac{2\pi^2}{3}-6
\end{align}
Therefore,
\begin{align}
\int^1_0\ln^2{x}\ln^2(1-x)\ {\rm d}x
&=2\left[-\frac{\pi^4}{30}-2\zeta(3)+\frac{91}{16}+\frac{5}{16}+\frac{\pi^4}{36}-2\zeta(3)-\frac{2\pi^2}{3}+6\right]\\
&=-\frac{\pi^4}{90}-8\zeta(3)-\frac{4\pi^2}{3}+24
\end{align}
A: Starting with the Beta function
\begin{align}
B(x, y) = \int_{0}^{1} t^{x-1} (1-t)^{y-1} \, dt
\end{align}
differentiate with respect to $x$ and $y$ twice. This leads to
\begin{align}
I &= \int_{0}^{1} t^{x-1} (1-t)^{y-1} \, \ln^{2}(t) \, \ln^{2}(1-t) \, dt \\
&= \partial_{x}^{2} \partial_{y}^{2}  B(x,y).
\end{align} 
Now the integral in question is the case for $x=y=1$. For this, it is seen that
\begin{align}
\partial_{x}^{2} \left. B(x,y) \right|_{x=1} = \frac{1}{y} \left[ \psi'(1) - \psi'(y+1) + \left( \psi(1) - \psi(y+1) \right)^{2} \right].
\end{align}
Differential with respect to $y$ yields
\begin{align}
\partial_{y}^{2} \partial_{x}^{2} \left. B(x,y) \right|_{x=1} &= \frac{1}{y} \left[ - \psi'''(y+1) - 2 \psi''(y+1) \left( \psi(1) - \psi(y+1) \right) + 2 \left( \psi'(y+1) \right)^{2} \right] \\
& \hspace{5mm} - \frac{2}{y^{2}} \left[ - \psi''(y+1) - 2 \psi'(y+1) \left( \psi(1) - \psi(y+1) \right) \right] \\
& \hspace{10mm} + \frac{2}{y^{3}} \left[ \psi'(1) - \psi'(y+1) + \left( \psi(1) - \psi(y+1) \right)^{2} \right].
\end{align}
Now
\begin{align}
\partial_{y}^{2} \partial_{x}^{2} \left. B(x,y) \right|_{x=y=1} &= 4 - \psi'''(2) + 4 \psi''(2) - 4 \psi'(2) + 2 \left( \psi'(2) \right)^{2} .
\end{align}
Since 
\begin{align}
\psi'(2) &= \frac{\pi^{2}}{6} -1  \\
\psi''(2) &= 2 - 2\zeta(3) \\
\psi'''(2) &= \frac{\pi^{4}}{45} - 6
\end{align}
then
\begin{align}
\partial_{y}^{2} \partial_{x}^{2} \left. B(x,y) \right|_{x=y=1} = 24 - \frac{4 \pi^{2}}{3} - \frac{\pi^{4}}{90} - 8 \zeta(3).
\end{align}
