# Expected number of lines meeting four given lines or “what is 1.72…”

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?!

• Note that the integrand is periodic in each direction with period pi, so you can cut your amount of work by a factor of 64. There's another factor of 6 from the symmetry of permuting (s1, t1), (s2, t2) and (s3, t3). A preliminary run using Gauss-Legendre quadrature with ~ 99^6/6 points gave an approximation of 1.726209603407 within 12 seconds; I'll try a larger net (maybe 199^6) reasonably soon. – Adam P. Goucher Jan 26 '17 at 19:47
• 1.726228333055 using 199^6 points. Now trying 299^6 points... – Adam P. Goucher Jan 26 '17 at 20:07
• @AdamP.Goucher: I didn't use a grid, but choose the points uniformly random, so I guess it doesn't make much of a difference?! – Moritz Firsching Jan 26 '17 at 21:51
• I'm using Gauss-Legendre quadrature, so the grid is non-uniform. It's unlikely (but not impossible) to be any worse than Monte Carlo, and can sometimes be significantly better (e.g. for smooth functions). Also, this works out faster since I don't need to generate random numbers or perform trig functions on the fly. By the way, my latest result is 1.726230356553 with 299^6 points (in 95 minutes). I'm doing an overnight run with 401^6 points, starting now... – Adam P. Goucher Jan 26 '17 at 21:58
• @AdamP.Goucher: great, thanks! that should give a few more digits.. – Moritz Firsching Jan 26 '17 at 22:03

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.

• A perspicuous algebraic form for the integral is: $$2^{-13}\int_{[-\infty,\infty]^6} \frac{\big|uv-vw+wx-xy+yz-zu \big|}{\big( (1+u^2)(1+v^2)(1+w^2)(1+x^2)(1+y^2)(1+z^2) \big)^{3/2}}du\,dv\,dw\,dx\,dy\,dz$$ which uses $u=\cot(t1),\ v=\cot(s2),\ w=\cot(t3),\ x=\cot(s1),\ y=\cot(t2),\ z=\cot(s3)$ – Matt F. Jan 29 '17 at 17:24
• nice work, Adam! @MattF.: How did you find this algebraic form? – Moritz Firsching Jan 30 '17 at 12:03
• @MoritzFirsching: first factor out two sins from each column, leaving only 1's and cot's in the determinant; then write sin as 1/sqrt(1+cot^2); and then multiply by some more of those expressions to account for dt1/du, etc. Finally make the coefficient 2^-7, since u runs through the real line twice as t1 goes around the circle. – Matt F. Jan 30 '17 at 17:54
• @MattF. Presumably the $2^{-13}$ is meant to be $2^{-7}$, then? – Adam P. Goucher Jan 31 '17 at 0:23
• @AdamP.Goucher That would be great! – Moritz Firsching Jan 31 '17 at 9:10