Solutions for 16 October 2026

Write \(I\) for the integral. We show that \(I=\tfrac{\pi^2}{8}\approx1.2337006\).

Why the obvious routes stall. Integrating in \(k\) first is elementary, but for each \(\theta\) it produces a mixture of arctangents and logarithms of \(\sin\theta\) and \(\cos\theta\) whose \(\theta\)-integral is unpleasant. Integrating in \(\theta\) first gives \(K(k)=\int_0^{\pi/2}\frac{d\theta}{\sqrt{1-k^2\sin^2\theta}}\), the complete elliptic integral of the first kind, which has no elementary antiderivative in \(k\). So \(I=\int_0^1\frac{K(k)}{1+k}\,dk\), and the weight \(\frac{1}{1+k}\) is the signal: it is exactly the factor that appears in Landen's transformation.

Step 1: Landen's transformation. For \(0\le k<1\) put \(\ell=\frac{2\sqrt k}{1+k}\). We claim \(K(\ell)=(1+k)K(k)\). By Gauss, \(K(k)=\frac{\pi}{2\,M(1,k')}\), where \(k'=\sqrt{1-k^2}\) and \(M\) is the arithmetic-geometric mean. Here \(1-\ell^2=\bigl(\frac{1-k}{1+k}\bigr)^2\), so \(\ell'=\frac{1-k}{1+k}\), and by homogeneity \(M(1,\ell')=\frac{1}{1+k}M(1+k,1-k)\). One AGM step sends \((1+k,1-k)\) to \(\bigl(1,\sqrt{1-k^2}\bigr)\), so \(M(1+k,1-k)=M(1,k')\). Hence \(K(\ell)=\frac{\pi(1+k)}{2M(1,k')}=(1+k)K(k)\).

Step 2: change the modulus. By Step 1, \(I=\int_0^1\frac{K(\ell)}{(1+k)^2}\,dk\). As \(k\) runs from \(0\) to \(1\), \(\ell\) increases from \(0\) to \(1\), and \(\frac{d\ell}{dk}=\frac{1-k}{\sqrt k\,(1+k)^2}\), so \(\frac{dk}{(1+k)^2}=\frac{\sqrt k}{1-k}\,d\ell\). Since \(\frac{\ell}{\sqrt{1-\ell^2}}=\frac{2\sqrt k}{1-k}\), this is \(\frac12\cdot\frac{\ell\,d\ell}{\sqrt{1-\ell^2}}\). Therefore

Step 3: back to a double integral, the other way round. Write \(K(\ell)\) as its \(\theta\)-integral and swap the order (the integrand is positive, so Tonelli applies): \(I=\frac12\int_0^{\pi/2}J(\sin\theta)\,d\theta\), where \(J(s)=\int_0^1\frac{\ell\,d\ell}{\sqrt{1-\ell^2}\sqrt{1-s^2\ell^2}}\). Put \(v=\ell^2\) and then \(w=\sqrt{1-v}\): \(J(s)=\int_0^1\frac{dw}{\sqrt{1-s^2+s^2w^2}}=\frac1s\operatorname{arsinh}\frac{s}{\sqrt{1-s^2}}\). Since \(\operatorname{arsinh}\frac{s}{c}=\ln\frac{1+s}{c}\) with \(c=\sqrt{1-s^2}\), this equals \(\frac{1}{2s}\ln\frac{1+s}{1-s}=\frac{\operatorname{artanh}s}{s}\) for \(0<s<1\).

Step 4: Wallis and the arcsine series. Expand \(\frac{\operatorname{artanh}s}{s}=\sum_{n\ge0}\frac{s^{2n}}{2n+1}\). All terms are positive, so we may integrate term by term, and Wallis gives \(\int_0^{\pi/2}\sin^{2n}\theta\,d\theta=\frac{\pi}{2}\binom{2n}{n}4^{-n}\). Hence

The series is \(\arcsin x=\sum_{n\ge0}\binom{2n}{n}\frac{4^{-n}x^{2n+1}}{2n+1}\) at \(x=1\). Its terms are of order \(n^{-3/2}\), so it converges at \(x=1\), and by Abel's theorem its sum is \(\arcsin1=\frac{\pi}{2}\). Therefore \(I=\frac{\pi}{4}\cdot\frac{\pi}{2}=\frac{\pi^2}{8}\).

Check. We evaluated the original double integral by two-dimensional composite Gauss-Legendre quadrature (100 nodes, 40 panels in \(\theta\), inner variable \(k=1-v^4\) to tame the endpoint) and agreed with \(\pi^2/8\) to \(4\times10^{-16}\). The single integral \(\int_0^1\frac{K(k)}{1+k}dk\), with \(K\) computed by the AGM, agreed to \(2\times10^{-16}\) on 40 and 200 panels.

The answer is \(\frac{e^{-2\gamma}}{2}\approx0.157618\), where \(\gamma\) is Euler's constant.

Why the obvious route stalls. One can write \(P(r_1+\dots+r_n<1)\) as an \(n\)-fold integral and let \(n\to\infty\), but these integrals show no usable pattern. The distribution of the full sum is in fact piecewise transcendental beyond \(1\). The way through is to notice that on \([0,1]\) its distribution function has a very simple shape with a single unknown constant, and then to find that constant from the behaviour of the Laplace transform at infinity.

Step 0: the law of one step. Given \(r_{k-1}\), the chance that the new point lies within distance \(xr_{k-1}\) of the origin is the area ratio \(x^2\) for \(0\le x\le1\). So \(r_k=r_{k-1}X_k\), where \(X_1,X_2,\dots\) are independent with density \(2x\) on \([0,1]\), and \(r_k=X_1X_2\cdots X_k\). Since \(E\,r_k=(2/3)^k\), the sum \(S=\sum_{k\ge1}r_k\) has \(E\,S=2\) and is finite almost surely. Moreover \(S=X_1(1+S')\), where \(S'=\sum_{k\ge2}X_2\cdots X_k\) has the same law as \(S\) and is independent of \(X_1\).

Step 1: the shape on \([0,1]\). Let \(F(t)=P(S<t)\). For \(0\le t\le1\) we have \(\frac{t}{1+S'}\le1\), so conditioning on \(S'\) gives \(F(t)=E\bigl[P\bigl(X_1<\tfrac{t}{1+S'}\bigm|S'\bigr)\bigr]=E\bigl[\tfrac{t^2}{(1+S')^2}\bigr]=ct^2\), with \(c=E\bigl[(1+S)^{-2}\bigr]\). The required probability is \(F(1)=c\).

Step 2: the Laplace transform. Let \(\varphi(\lambda)=E\,e^{-\lambda S}\) for \(\lambda\ge0\). It is continuous, positive, and \(\varphi(0)=1\). From \(S=X_1(1+S')\),

The right side is continuously differentiable, hence so is \(\varphi\) on \((0,\infty)\), and differentiating gives \(2\lambda\varphi+\lambda^2\varphi'=2\lambda e^{-\lambda}\varphi\), that is \(\frac{\varphi'}{\varphi}=-\frac{2(1-e^{-\lambda})}{\lambda}\). The right side is bounded near \(0\), so integrating from \(0\) gives

Step 3: the behaviour at infinity. Integration by parts gives \(\int_0^\lambda\frac{1-e^{-u}}{u}\,du=(1-e^{-\lambda})\ln\lambda-\int_0^\lambda e^{-u}\ln u\,du\), with no boundary term at \(0\). As \(\lambda\to\infty\) the last integral tends to \(\Gamma'(1)=-\gamma\), so the left side equals \(\ln\lambda+\gamma+o(1)\). Hence \(\lambda^2\varphi(\lambda)\to e^{-2\gamma}\).

Step 4: the same limit from Step 1. Split the transform at \(t=1\): \(\varphi(\lambda)=\int_0^1 2ct\,e^{-\lambda t}\,dt+\int_{[1,\infty)}e^{-\lambda t}\,dF(t)\). The first term is \(\frac{2c}{\lambda^2}\bigl(1-(1+\lambda)e^{-\lambda}\bigr)\) and the second lies between \(0\) and \(e^{-\lambda}\). Hence \(\lambda^2\varphi(\lambda)\to2c\). Comparing with Step 3, \(2c=e^{-2\gamma}\), so the probability is \(c=\frac{e^{-2\gamma}}{2}\).

Remark. \(S\) is a generalised Dickman variable with parameter \(2\). Choosing the point uniformly on a segment instead of a disc gives the classical Dickman law, and the same argument then yields \(e^{-\gamma}\).

Check. A Monte Carlo run of \(4\times10^7\) sequences (80 terms each) gave \(0.157539\), within \(1.4\) standard errors of \(0.157618\), and the ratio \(P(S<x)/x^2\) was flat on \([0,1]\). Integrating the delay equation \(xf'(x)=f(x)-2f(x-1)\) for the density out to \(x=14\) reproduced the constant \(e^{-2\gamma}\) to \(10^{-9}\).