Expected number of lines meeting four given lines or "what is 1.72..."

The integrand is periodic modulo $\pi$ in each variable, so it suffices to integrate each variable over $[0, \pi]$ and replace the constant factor by $2^{-7}$.

If we were to apply a change of variables (e.g. set $x = \cos(s_1)$ and similarly for the other five variables), we would have an integral of a piecewise-algebraic function which thus belongs to Kontsevich and Zagier's 'ring of periods'. Alas, there is no proven algorithm for determining whether periods are expressible in terms of elementary functions.

Now, using Gauss-Legendre integration with $N$ points in each variable, we can approximate it with a weighted sum of $N^6$ evaluations of the integrand. Since the absolute value of the determinant is unchanged under permuting the columns, and zero whenever two columns are equal, we can reduce this to $\binom{N^2}{3}$ function evaluations.

I managed to implement this in C code, where each iteration, amortized, only takes four floating-point additions, three multiplications, and an absolute-value calculation. With several hours running on a 144-core machine, it could do the calculation for both $N = 401$ and $N = 409$. The results were 1.726230867 and 1.726230885, respectively.

Consequently, the integral is roughly 1.726230876, with an expected error on the order of $10^{-8}$.

---

EDIT: Matt F. posted an algebraic form for the integral, by taking $u = \cot(t1)$, $v = \cot(s2)$, $w = \cot(t3)$, $x = \cot(s1)$, $y = \cot(t2)$, $z = \cot(s3)$:

$$ I = 2^{-7} \int \dfrac{|uv-vw+wx-xy+yz-zu|}{[(1+u^2)(1+v^2)(1+w^2)(1+x^2)(1+y^2)(1+z^2)]^{3/2}} $$

where the integral is taken over $(u,v,w,x,y,z) \in \mathbb{R}^6$. Observing that the numerator is twice the area of a triangle with vertices $(u, x), (y, v), (w, z)$, we can perform a change of variables by setting:

$$ w = u + (y-u)c - (v-x)h \textrm{ ; } z = x + (v-x)c + (y-u)h $$

where $c \in \mathbb{R}$ and $h \in \mathbb{R}^{+}$. Taking the Jacobian into account, the integral becomes:

$$ I = 2^{-6} \int \dfrac{h[(y-u)^2+(v-x)^2]^2}{[(1+u^2)(1+y^2)(1+(u + (y-u)c - (v-x)h)^2)(1+x^2)(1+v^2)(1+(x + (v-x)c + (y-u)h)^2)]^{3/2}} $$

where $(u,y,x,v,c,h) \in \mathbb{R}^5 \times \mathbb{R}^+$. Although the integral is uglier, it satisfyingly has no absolute value operator.

If anyone is still interested, Antonio Lerario and I recently published a paper: Probabilistic Schubert Calculus: Asymptotics, in which we give a more convenient formula to compute this number.

What was called before $\operatorname{edeg}G(2,4)$ is now denoted by $\delta_{1,3}$ and we give a line integral formula in Proposition 24:

\begin{aligned} \delta _{1,3}=-6 \pi ^{4}\int _{0}^1 L(u)^{2}\mathrm {sinh}(w(u))w'(u)\mathrm{d}u \end{aligned} where $L=F⋅G$ and $w=log(F/G)$ with \begin{aligned} F(u)&:=\int _{0}^{\pi /2} \frac{u \ \sin ^2(\theta )}{\sqrt{\cos ^2\theta +u^2 \sin ^2\theta }}\mathrm {d}\theta \\ G (u)&:=\int _{0}^{\pi /2} \frac{\sin ^2(\theta )}{\sqrt{\sin ^2\theta +u^2 \cos ^2\theta }}\mathrm {d}\theta . \end{aligned}

We didn't run any advanced numerical evaluation but this looks much nicer than the previous formula. I hope this will help. It is still unknown to us if this number can be expressed as a closed formula using special functions.