Calculate in closed form $\int_0^1 \int_0^1 \frac{dx\,dy}{1-xy(1-x)(1-y)}$ Can we possibly compute the following integral in terms of known constants?
$$\int_0^1 \int_0^1 \frac{dx\,dy}{1-xy(1-x)(1-y)}$$
Some progress was already done here http://integralsandseries.prophpbb.com/topic279.html but still we have a hypergeometric function. What's your thoughts on it? 
UPDATE: The question was also posted on mathoverflow here https://mathoverflow.net/questions/206253/calculate-in-closed-form-int-01-int-01-fracdx-dy1-xy1-x1-y#
 A: Let $x=\sin^2(t)$ and $y=\sin^2(w)$, we obtain
\begin{align}
I & = \int_0^{\pi/2}\int_0^{\pi/2} \dfrac{4\sin(t)\cos(t)\sin(w)\cos(w)dtdw}{1-\sin^2(t)\cos^2(t)\sin^2(w)\cos^2(w)}\\
& = 4\sum_{k=0}^{\infty}\int_0^{\pi/2}\int_0^{\pi/2} \sin^{(2k+1)}(t) \cos^{(2k+1)}(t)\sin^{(2k+1)}(w) \cos^{(2k+1)}(w)dtdw\\
& = 4 \sum_{k=0}^{\infty} \left(\int_0^{\pi/2}\sin^{(2k+1)}(t) \cos^{(2k+1)}(t)dt\right)^2\\
& = 4 \sum_{k=0}^{\infty} \dfrac1{4^{2k+1}}\left(\int_0^{\pi/2}\sin^{(2k+1)}(2t)dt\right)^2 = \sum_{k=0}^{\infty} \dfrac1{4^{2k}}\left(\int_0^{\pi}\sin^{(2k+1)}(t)\dfrac{dt}2\right)^2\\
& = \sum_{k=0}^{\infty} \dfrac1{4^{2k}}\left(\int_0^{\pi/2}\sin^{(2k+1)}(t)dt\right)^2
\end{align}
Recall that $$\int_0^{\pi/2}\sin^{(2k+1)}(t)dt=\dfrac{4^k}{2k+1}\dfrac1{\dbinom{2k}k}$$
This gives us
$$I = \sum_{k=0}^{\infty} \dfrac1{(2k+1)^2} \dfrac1{\dbinom{2k}k^2} = _3F_2\left(1,1,1;\frac32,\frac32;\frac1{16}\right)$$
A: An incomplete answer.
Tackling the inner integral with elementary methods,
\begin{align}
\int^1_0\frac{{\rm d}x}{1-ax+ax^2}
&=\int^1_0\frac{{\rm d}x}{a\left(x-\frac{1}{2}\right)^2+1-\frac{a}{4}}\\&=\int^\frac{1}{2}_{-\frac{1}{2}}\frac{{\rm d}x}{ax^2+1-\frac{a}{4}}\\
&=\frac{1}{\sqrt{a\left(1-\frac{a}{4}\right)}}\int^\frac{\sqrt{a}}{2\sqrt{(1-\frac{a}{4})}}_{-\frac{\sqrt{a}}{2\sqrt{(1-\frac{a}{4})}}}\frac{{\rm d}x}{1+x^2}\\&=\frac{4}{\sqrt{a(4-a)}}\arctan\left(\sqrt{\frac{a}{4-a}}\right)
\end{align}
Thus
\begin{align}
\iint\limits_{[0,1]^2}\frac{{\rm d}A}{1-xy(1-x)(1-y)}
&=4\int^1_0\arctan\left(\sqrt{\frac{x(1-x)}{x^2-x+4}}\right)\frac{{\rm d}x}{\sqrt{x(1-x)(x^2-x+4)}}\\
&=16\int^1_0\arctan\left(\sqrt{\frac{1-x^2}{15+x^2}}\right)\frac{{\rm d}x}{\sqrt{(1-x^2)(15+x^2)}}\\
&=\frac{16}{\sqrt{15}}\int^1_0\operatorname{F}\left(\arcsin{x}\right)\frac{x\ {\rm d}x}{\sqrt{(1-x^2)(15+x^2)}}\\
&=\frac{8}{15}\int^{\operatorname{K}}_{-\operatorname{K}}u\operatorname{sn}(u)\ {\rm d}u
\end{align}
where $\operatorname{sn}(u)$ is a Jacobi elliptic function with modulus $k=15^{-1/2}i$. 
Let $R$ be the parallelogram with vertices $-\operatorname{K}-2i\operatorname{K}'$, $\operatorname{K}-2i\operatorname{K}'$, $\operatorname{K}+2i\operatorname{K}'$, $-\operatorname{K}+2i\operatorname{K}'$. Parameterising along the contour,
\begin{align}
\int_{\partial R}z^{\color{red}{2}}\operatorname{sn}(z)\ {\rm d}z
&=-8i\operatorname{K}\left(\frac{4}{\sqrt{15}}\right)\int^{\operatorname{K}}_{-\operatorname{K}}u\operatorname{sn}(u)\ {\rm d}u+2i\int^{2\operatorname{K}'}_{-2\operatorname{K}'}(\operatorname{K}^2-t^2)\operatorname{cd}(it)\ {\rm d}t\\
&=2\pi i\sum_{\pm}\operatorname*{Res}_{z=\pm i\operatorname{K}'}z^2\operatorname{sn}(z)=-4\pi\sqrt{15}\operatorname{K}\left(\frac{4}{\sqrt{15}}\right)^2
\end{align}
$\operatorname{cd}(t)$ has an anti-derivative, namely
$$2i\operatorname{K}^2\int^{2\operatorname{K}'}_{-2\operatorname{K}'}\operatorname{cd}(it)\ {\rm d}t=-4i\sqrt{15}\operatorname{K}^2\left.\ln\left(\operatorname{nd}(t)+\frac{i}{\sqrt{15}}\operatorname{sd}(t)\right)\right|^{2i\operatorname{K}'}_0$$
Otherwise, this seems to be as far as I can go for now.
A: From hypergeometric function to elliptic integrals
The hypergeometric function $_3F_2$ from your link (or user17762's answer) may be rewritten (using following integral representation and $_2F_1\left(1,1;\frac32; t\right)=\frac d{dt}\left[\left(\arcsin\sqrt{t}\,\right)^2\right]$) as :
\begin{align}
\tag{1}I&=\;_3F_2\left(1,1,1;\frac32,\frac32;\frac1{16}\right)=\frac 12\int_0^1\frac{_2F_1\left(1,1;\frac32;\frac x{16}\right)}{\sqrt{1-x}}dx\\
&=2\int_0^1\frac{\arcsin\frac{\sqrt{x}}4}{\sqrt{x\,(1-x)\left(1-\frac x{16}\right)}}dx\\
\tag{2}&=16\int_0^{1/4}\frac{\arcsin u}{\sqrt{(1-16 u^2)\left(1-u^2\right)}}\,du\\
&\tag{3}=16\int_0^{\arcsin(1/4)} \frac {y}{\sqrt{1-4^2\sin(y)^2)}}\;dy\\
&=\left.16\;y\;F\left(y, 4\right)\right|_0^{\arcsin(1/4)}-16\int_0^{\arcsin(1/4)}\;F\left(y, 4\right)\;dy\\
\tag{4}&=4\;\arcsin(1/4)\;K\left(\frac 1{4}\right)-16\int_0^{\arcsin(1/4)}\;F\left(y, 4\right)\;dy\\
\end{align}
With $\;\displaystyle F(\phi,k):=\int_0^{\phi} \frac {dy}{\sqrt{1-k^2\sin(y)^2)}},\;K(k):=F\left(\frac{\pi}2,k\right)\;$ the incomplete and complete elliptic integral of the first kind $F$ and $K$ (one should take care of replacing $F(\phi,k)$ by EllipticF[$\phi,k^2$] and $K(k)$ by EllipticK[$k^2$] while using Alpha/Mathematica).
'$I$' may be expressed with elliptic integrals but this doesn't seem to help for a closed form... Observing that $\;\displaystyle K\left(\frac 1{4}\right)=\frac{\pi}2\sum_{k=0}^\infty\frac{\binom{2k}{k}^2}{2^{8\,k}}\;$ let's start it all again using the 
Squared central binomial series
We want $\;\displaystyle I=\sum_{k=0}^\infty \frac 1{(2k+1)^2\binom{2k}{k}^2}\;$ but Nected Batir proved in 2004 following more general formula (for $n\in\mathbb{N},\; n\ge3$) :
\begin{align}
\tag{5}I_n(x)&:=\sum_{k=0}^\infty \frac{x^{2k+1}}{(2k+1)^n\,\binom{2k}{k}^2}\\
&=\frac{4\,(-1)^{n-3}}{(n-3)!}\int_0^{\pi/2}\int_0^{\arcsin(x\cos(y)/4)}\frac t{\sin(t)}\;\log^{n-3}\left(\frac{4\sin(t)}{x\cos(y)}\right)\;dt\,dy\\
\end{align}
For $\,n=3\,$ this is simply
$\;\displaystyle I_3(x)=\int_0^{\pi/2}\int_0^{\arcsin(x\cos(y)/4)}\frac t{\sin(t)}\;dt\,dy\;\;$  and after derivation
\begin{align}
I_2(x)&=x\,I_3(x)'=4\,x\int_0^{\pi/2} \frac {\arcsin(x\cos(y)/4)}{x\cos(y)/4} \frac {\partial}{\partial x}\left[\arcsin(x\cos(y)/4)\right]\;dy\\
\tag{6}I_2(x)&=4\int_0^{\pi/2} \frac {\arcsin(x\cos(y)/4)}{\sqrt{1-(x\cos(y)/4)^2}} \;dy\\
\end{align}
Unfortunately this is merely $(2)$ in the case $x=1$ (the integral for $x=4$ is divergent and $I_3(4)=8\pi G-14\zeta(3)\;$ but this won't help here...) so nothing complete here either...
The journey may be of interest anyway...
A: Our integral equals:
$$\begin{eqnarray*} I=\sum_{n\geq 0}\iint_{[0,1]^2}x^n y^n(1-x)^n (1-y)^n\,dx\,dy &=& \sum_{n\geq 0}\left(\int_{0}^{1}x^n(1-x)^n\,dx\right)^2\\&=&\sum_{n\geq 0}\left(\frac{\Gamma(n+1)^2}{\Gamma(2n+2)}\right)^2\\&=&\sum_{n\geq 0}\frac{1}{(2n+1)^2\binom{2n}{n}^2}\tag{1}\end{eqnarray*}$$
and:
$$ \sum_{n\geq 0}\frac{x^n}{(2n+1)\binom{2n}{n}}=\frac{4\arcsin\frac{\sqrt{x}}{2}}{\sqrt{x(4-x)}}\tag{2}$$
hence $I$ is the squared $L^2$-norm of a quite complicated function:
$$ I = \frac{16}{\pi}\int_{0}^{\pi}\frac{\arcsin\frac{e^{it/2}}{2}\cdot\arcsin\frac{e^{-it/2}}{2}}{\sqrt{(4-e^{it})(4-e^{-it})}}\,dt.\tag{3}$$
