Sylvester’s four-point problem in a disc, Geometric Probability

Problem 46.1

Let $P_1, P_2, P_3, P_4$ be four points chosen independently and uniformly at random in a disc of radius $r$ . What is the probability that they form the vertices of a convex quadrilateral?

Solution. By the same dichotomy as in Chapters 43 and 45, $\mathbb{P}(\text{convex}) = 1 - 4 \cdot \frac{\mathbb{E}[\operatorname{Area}(\triangle P_1 P_2 P_3)]}{\pi r^2}.$

The expected triangle area for three uniform random points in a disc of radius $r$ is the classical constant $\mathbb{E}[\operatorname{Area}(\triangle P_1 P_2 P_3)] = \frac{35}{48 \pi} \cdot r^2.$ This is due to Woolhouse (1867) and can be derived by polar-coordinate integration: conditioning on the distances of two of the three points from the origin, one reduces the calculation to a one-dimensional integral that eventually yields the expected area $\tfrac{35}{48\pi}\,r^2$ . As a fraction of the disc’s own area $\pi r^2$ , this is $\tfrac{35}{48\pi^2} \approx 0.0739$ (not $\tfrac{35}{48\pi}$ , which is the area itself when $r = 1$ ). (The computation is longer than those for the triangle and square cases, but entirely elementary.)

Substituting, $\mathbb{P}(\text{convex}) = 1 - 4 \cdot \frac{(35/(48\pi)) \cdot r^2}{\pi r^2} = 1 - \frac{35}{12 \pi^2}.$ Hence $\boxed{\mathbb{P}(\text{convex}) = 1 - \frac{35}{12\pi^2} \approx 0.7045.} \qedhere$

Compared with the triangle ( $\tfrac23$ ) and the square ( $\tfrac{25}{36}$ ), the disc gives the largest convex-position probability among the three. Blaschke’s 1917 theorem says that among all convex regions of area $1$ , the triangle minimises this probability and the ellipse (or disc) maximises it. The extremal values $\tfrac23$ and $1 - \tfrac{35}{12\pi^2}$ frame the entire range of Sylvester’s four-point problem over convex shapes.

import numpy as np
N = 10**6
# Uniform in unit disc via (sqrt(u), angle)
r = np.sqrt(np.random.rand(4, N))
t = np.random.rand(4, N) * 2*np.pi
u, v = r*np.cos(t), r*np.sin(t)
def inside_tri(pu, pv, au, av, bu, bv, cu, cv):
    det = (bv-cv)*(au-cu) + (cu-bu)*(av-cv)
    w1 = ((bv-cv)*(pu-cu) + (cu-bu)*(pv-cv)) / det
    w2 = ((cv-av)*(pu-cu) + (au-cu)*(pv-cv)) / det
    return (w1 > 0) & (w2 > 0) & (1 - w1 - w2 > 0)
any_in = np.zeros(N, dtype=bool)
for i in range(4):
    j, k, l = [m for m in range(4) if m != i]
    any_in |= inside_tri(u[i], v[i], u[j], v[j], u[k], v[k], u[l], v[l])
print(f"sim: {(1 - any_in).mean():.5f}   exact: {1 - 35/(12*np.pi**2):.5f}")
# sim: 0.70427   exact: 0.70448