Given four random lines in $\mathbb{R}P^3$, how many lines intersect all of those lines?
In the recent paper [Probabilistic Schubert Calculus](http://arxiv.org/abs/1612.06893v1), Peter Bürgisser and Antonio Lerario
discuss this question and much more general versions of it.

In **Proposition 6.7**, they determine the expected number of lines meeting four random lines in $\mathbb{R}P^3$ to be:

 $$\operatorname{edeg}G(2,4) =\\ 2^{-13}\int_{[0,2\pi]^6}\left|
 \det\begin{pmatrix}
 \sin{t_1}\sin{s_1}&\sin{t_2}\sin{s_2}&\sin{t_3}\sin{s_3}\\
 \cos{t_1}\sin{s_1}&\cos{t_2}\sin{s_2}&\cos{t_3}\sin{s_3}\\
 \sin{t_1}\cos{s_1}&\sin{t_2}\cos{s_2}&\sin{t_3}\cos{s_3}
 \end{pmatrix}\right|dt_1dt_2dt_3ds_1ds_2ds_3$$ 

The integrand can be expanded to
$$\left|\cos(s_2)\sin(s_1)\sin(s_3)\sin(t_2)\sin(t_1 - t_3)- \sin(s_2)\big(\cos(s_1)\sin(s_3)\sin(t_1)\sin(t_2 - t_3) + \cos(s_3)
\sin(s_1)\sin(t_3)\sin(t_1 - t_2)\big)\right|$$

The absolute value in the integrand and the integration over six variables seem to make it difficult to evaluate this integral to high precision.
I am interested in the exact value of this number. Bürgisser and Lerario computed it to be $1.72$ rounded to two digits;
I ran $10^{11}$ Monte Carlo evaluations to obtain $1.726225\dots$
with an estimated error of $7.3\cdot 10^{-6}$

My questions are:

 - Can the exact value of $\operatorname{edeg}G(2,4)$ be determined; is there perhaps a relation with other constants; is it, for example, algebraic over $\mathbb{Q}[\pi]$?
 - If this is too difficult, what are ways of evaluating the above integral to higher accuracy numerically?

**Update**: In AdamP.Goucher's great answer, he provides more digits ($1.726230876$) and also uses a reformulation described by Matt F. in a comment to get an integrand without any absolute values:

$$G(2,4) = 2^{-6} \int \dfrac{h\left((y-u)^2+(v-x)^2\right)^2}{\left((1+(u + (y-u)c - (v-x)h)^2)(1+(x + (v-x)c + (y-u)h)^2)\prod_{\alpha\in\{u,y,x,v\}}(1+\alpha^2)\right)^{3/2}}$$
where we integrate over $\mathbb{R}^5\times \mathbb{R}^+$ (and $h$ is the positive variable).

Perhaps this form could come in handy for evaluating this integral to higher accuracy?!

--------

**Update 2** The new formula from L.Mathis and Antonio Lerario is very useful for calculating digits! The following [mpmath][1] code can returns in less than $3$ minutes

 $1.7262312489219034885256331685361697650475579915479447$ 

(the last few digits might not be accurate)
I expect to make that even much faster when solving the two integrals $F$ and $G$ symbolically first. 

    import functools
    from mpmath import mp
    @functools.lru_cache(maxsize=1000)
    def F(u):
        return mp.quad(lambda phi: (u*mp.sin(phi)**2)/(mp.sqrt(mp.cos(phi)**2 + u**2*mp.sin(phi)**2)), [0, mp.pi/2])
    @functools.lru_cache(maxsize=1000)
    def G(u):
        return mp.quad(lambda phi: (mp.sin(phi)**2)/(mp.sqrt(mp.sin(phi)**2 + u**2*mp.cos(phi)**2)), [0, mp.pi/2])
    def H(u):
        return F(u)/G(u)
    def L(u):
        return F(u)*G(u)
    def integrand3(u):
        return L(u)**2*(1/H(u) - H(u))*mp.diff(H, u)/H(u)

    dps = 50
    mp.dps = dps

    %time z = 3*mp.quad(integrand, [0,1]); z


  [1]: http://mpmath.org/