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

Question

Given four random lines in $\mathbb{R}P^3$, how many lines intersect all of those lines? In the recent paper Probabilistic Schubert Calculus, 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 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

Answer 1

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.

Answer 2

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.