Aquí está mi computación. En primer lugar, no tenía ni idea de cómo integrar funciones trigonométricas implícitas, así que decidí cambiar a la integración con respecto a los lados $x,y,z$ del triángulo $XYZ$. Recuerde que las condiciones para el dominio son $xz\le y^2+xy$, $x>y$ más la desigualdad del triángulo $z>x-y$. Esta desigualdad del triángulo es compatible con la primera condición si y solo si $x< (\sqrt 2+1)y$. Así que el dominio admisible está dado por $$ y< x< (\sqrt 2+1)y,\qquad x-y Ahora tenemos que integrar sobre los puntos en el círculo, es decir, para el caso cuando el radio circunscrito es $1$ y con respecto a las variables angulares. Relajamos el radio circunscrito a $r>0$, colocamos $Z$ en el semieje positivo y el centro del círculo en el origen (así $Z=r$), y tomamos $Y=re^{2i\varphi}, X=re^{2i\psi}$, entonces $$ x=2r\sin\varphi,\qquad y=2r\sin\psi,\qquad z=2r\sin(\varphi-\psi) $$ (el último hasta un signo) donde $\varphi,\psi$ se mueven independientemente sobre $[0,\pi]$. La condición $x>y$ reduce la medida $d\varphi\,d\psi$ a $\frac{\pi^2}2$ y cada terna admisible $x,y,z$ se realiza dos veces (hasta el volteo sobre los ejes horizontales). Podemos usar cualquier medida en $(0,+\infty)$ que queramos al integrar sobre $r$, así que elegiremos $\frac{dr}{r}$.
Ahora es hora de calcular el Jacobiano $J$ de la transformación $(r,\varphi,\psi)\mapsto (x,y,z)$. Tenemos $$ J=8\left|\begin{matrix} \sin\varphi & r\cos\varphi & 0 \\ \sin\psi & 0 & r\cos\psi \\ \sin(\varphi-\psi) &r\cos(\varphi-\psi) & -r\cos(\varphi-\psi) \end{matrix}\right| \\ =8r^2[\cos\varphi\cos\psi\sin(\varphi-\psi)-\sin\varphi\cos\psi\cos(\varphi-\psi)+\cos\varphi\sin\psi\cos(\varphi-\psi)] \\ =8r^2[\cos\varphi\cos\psi\sin(\varphi-\psi)- \sin(\varphi-\psi)\cos(\varphi-\psi)] \\ =8r^2\sin(\varphi-\psi)\sin\varphi\sin\psi=\frac{xyz}r\,, $$ es decir, $$ \frac{dr}r\,d\varphi\,d\psi=\frac{dx\,dy\,dz}{xyz}\,, $$ lo cual es bastante bueno.
Ahora permita que $\Omega$ sea el conjunto de todos los triángulos buenos y $\omega$ sea el conjunto de triángulos buenos con $r=1$. Entonces el cálculo de la servilleta es el siguiente: $$ \left[\int_0^\infty\frac{dr}{r}\right]\int_\omega d\varphi\,d\psi= 2\int_{\Omega}\frac{dr}{r}\,d\varphi\,d\psi =2\int_{\Omega}\frac{dx\,dy\,dz}{xyz} \\ = 2\int_0^\infty\frac{dy}{y}\left[\int_{y}^{(\sqrt 2+1)y}\frac{dx}x\left(\int_{x-y}^{y+\frac{y^2}x}\frac{dz}z\right)\right] $$ Sustituyendo $x=ty$, $t\in(1,\sqrt 2+1)$, evaluamos las integrales internas a $$ \int_1^{\sqrt 2+1}\log\frac{t+1}{t(t-1)}\frac{dt}t\,, $$ que es independiente de $y$, así que finalmente obtenemos $$2\left[\int_0^\infty\frac{dy}y\right] \int_1^{\sqrt 2+1}\log\frac{t+1}{t(t-1)}\frac{dt}t\,, $$ Ahora queda cancelar $\int_0^\infty\frac{ds}s$ para obtener la identidad reclamada.
El pequeño (pero algo irritante) problema es que la integral $\int_0^\infty \frac{ds}s$ diverge. Sin embargo, como solía decir Landau (el físico), "Un pollo no es un pájaro y un logaritmo no es infinito". Así que arreglaremos eso con un argumento un poco más elaborado. ${}{}$
Pondremos una condición adicional en los triángulos admisibles, a saber, exigiremos que $y$ sea comparable a $r$. Siempre trivialmente tenemos $r\ge \frac y2$, así que fijaremos un $A>0$ grande y exigiremos que $r. Note que es una condición solo sobre la forma, no sobre el tamaño. Denotando por $\omega_A$ el conjunto de triángulos admisibles con radio circunscrito $1$ y por $\Omega_A$ el conjunto de todos los triángulos admisibles y eligiendo un $R>0$ enorme, podemos escribir, como antes, $$ \left[\int_1^R\frac{dr}r\right]\int_{\omega_A}d\varphi\,d\psi=2\int_{\Omega_A\cap\{1 Ahora note que $$ \Omega_A\cap\{2 así que en lugar de una identidad, escribimos dos desigualdades y concluimos que $\log R\int_{\omega_A}d\varphi\,d\psi$ está atrapado entre $\log\frac{R}{2A}$ veces $2\int_1^{\sqrt 2+1}\left(\int_{E(t,A)}\frac{ds}{s}\right)\frac{dt}t$ y $\log(2AR)$ veces la misma cantidad, donde ponemos $x=ty, z=sy$ y definimos $E(t,A)$ como el conjunto de $s$ para el cual el triángulo con los lados $1,t,s$ es admisible.
Ahora, dividiendo por $\log R$, dejando que $R\to+\infty$ y usando el teorema del límite, concluimos que $$ \int_{\omega_A}d\varphi\,d\psi= 2\int_1^{\sqrt 2+1}\left(\int_{E(t,A)}\frac{ds}{s}\right)\frac{dt}t $$ para cada $A$ fijo. Pero cuando $A\to+\infty$, el conjunto $\omega_A$ se expande a $\omega$ y para cada $t$, el conjunto $E(t,A)$ se expande a $[t-1,1+\frac 1t]$, así que el teorema de convergencia monótona termina la historia.
Eso es. Intenté encontrar un cambio de variables simple que evitara usar identidades no triviales para $\mathrm{Li}_2$ para obtener la respuesta final de $\frac{\pi^2}8$, pero hasta ahora no lo logré. Sin embargo, tenemos el resultado, así que, espero, alguien eventualmente proporcionará una prueba más simple. El valor de $\frac 12$ de la probabilidad es demasiado agradable para estar allí sin una razón clara. :-)
