Here is the definition of the problem as phrased on Wikipedia:

Consider an equilateral triangle inscribed in a circle. Suppose a chord of the circle is chosen at random. What is the probability that the chord is longer than a side of the triangle?

The question is easy to visualize and might appear well-posed at first. For example, in the figure below, the green chord is longer that the side of the inscribed equilateral triangle, while the blue one is shorter. If one were to draw chords at random, what fraction of them would end up being longer than the side of the triangle?

In [4]:

```
R = 1 # radius of circle
```

In [5]:

```
# example plot
fig = plt.figure(figsize=(12,6)); ax = fig.add_subplot(1,1,1)
ax.set(xlim = (-2*R, 2*R), ylim = (-R - 0.1, R + 0.1), aspect = 1)
_tmp = ax.axis('off')
# circle
_tmp = ax.add_patch(patches.Circle((0,0), radius = R, fc = 'lightgrey', lw = 1))
# equilateral triangle inscribed to circle
_tmp = ax.add_patch(patches.RegularPolygon((0,0), 3, R, lw = 3, fill = False))
# add two chords, one shorter & one longer than triangle's side
_tmp = add_chord(ax, R, np.pi/3, -np.pi/12, "blue", 3)
_tmp = add_chord(ax, R, 5*np.pi/12, -2*np.pi/3, "green", 3)
```

*how exactly* a chord is chosen 'at random'. The 'paradox' in this problem is that different ways to choose a chord lead to different results. Bertrand considered three different methods to produce the random chords and obtained three different solutions.

*see E.T. Jaynes, "The well-posed problem", 1973*.

In [6]:

```
N = 100000 # number of simulations
triangle_edge_length = get_triangle_side_length(R) # side length of the inscribed equilateral triangle
```

The first method is the following.

"Choose two random points on the circumference of the circle and draw the chord joining them".

This is equivalent to selecting two random points $a$ and $b$ from the interval $[0, 2\pi)$. Given the center $c$ and radius $R$ of the circle, it is then straightforward to map $a$ and $b$ to two points on the circumference of the circle, calculate their cartesian coordinates and find the length of the chord they form.

For the purposes of simulation, we repeat the above process $N = 100,000$ times.

In [7]:

```
# each of the two endpoints of a single chord is
# defined by a (random) angle in [0, 2\pi)
end_a = 2 * np.pi * np.random.rand(N) # first
end_b = 2 * np.pi * np.random.rand(N)
# the cartesian (2d) coordinates for each endpoint
sol_A_coord_a = [(R*np.cos(theta), R*np.sin(theta)) for theta in end_a] # first endpoint
sol_A_coord_b = [(R*np.cos(theta), R*np.sin(theta)) for theta in end_b] # second endpoint
# the length of each chord
lengths = l2dist(sol_A_coord_a, sol_A_coord_b)
solution_A_pct = 100*sum(lengths > triangle_edge_length) / N
```

In [8]:

```
plot_length_distr_and_chords(lengths, R, solution_A_pct,
triangle_edge_length, sol_A_coord_a, sol_A_coord_b, N)
```

The second method proceeds differently.

"Choose a radius of the circle, choose a point on the radius and construct the chord through this point and perpendicular to the radius."

Selecting a radius of the circle is equivalent to selecting a random point $\phi \in [0, 2\pi)$. Moreover, selecting a random point on the radius is equivalent to selecting a random point $p \in [0, R)$. Given $\phi$ and $p$, it is straighforward to obtain the cartesian coordinates for the two endpoints $a$ and $b$ of the constructed chord.

In [9]:

```
phi = 2 * np.pi * np.random.rand(N) # choose random radius
point = R * np.random.rand(N) # choose random point on radius
# for the radius at \phi = 0, the two endpoints are positioned at
# angles +/- arccos(p/R)
absolute_angles = np.array([(np.arccos(p/R)) for p in point])
# selecting a radius at \phi, simply rotates the angles of
# endpoints by \phi, i.e. the two endpoints are
# positioned at phi +/- arccos(p/R)
sol_B_coord_a = [(R*np.cos(theta), R*np.sin(theta)) for theta in phi + absolute_angles]
sol_B_coord_b = [(R*np.cos(theta), R*np.sin(theta)) for theta in phi - absolute_angles]
# calculate the length of the constructed chord
lengths = l2dist(sol_B_coord_a, sol_B_coord_b)
solution_B_pct = 100*sum(lengths > triangle_edge_length)/N
```

In [10]:

```
plot_length_distr_and_chords(lengths, R, solution_B_pct,
triangle_edge_length, sol_B_coord_a, sol_B_coord_b, N)
```

We saw that the first method constructed chords by selecting random pairs of points from the circumference of the circle, while the second method selected random points from a radius of the circle. By contrast, this method proceeds by selecting random points from the *surface* within the circle (i.e., the chances that an area will incude one selected random point are proportional to its surface). Specifically, to follow this method,

"Choose a point anywhere within the circle and construct a chord with the chosen point as its midpoint."

Suppose we have selected a midpoint with the aforementioned method; let $r$ be its distance from the centre $c = (0,0)$ of the circle, and $s=s(r)$ be the surface enclosed in the 'small' circle with center $c$ and radius $r$.

Consider a newly drawn midpoint with distance $r'$ from center $c$. Observe that the probability that $r' < r$ is equal to $\pi r^2 / \pi R^2 = (r/R)^2$.

From the above observation, it is easy to see that selecting the midpoint is equivavalent to following these steps:

- select a random point $p \in [0,1)$,
- select the circle with center $c = (0,0)$ and radius $r = R * \sqrt{p}$,
- select a random point $m$ from its circumference (equivalently, select random $\phi \in [0, 2\pi)$,

- use $m$ as the method's midpoint.

Given $r$ and $\phi$ associated with the selected midpoint, it easy to calculate the cartesian coordinate of the two endpoints $a$ and $b$ of the constructed chord.

In [11]:

```
R_sq = np.square(R)
# pick midpoint so that
# midpoints are picked uniformly at random
# across the area inside the circle
mid_radius = R * np.sqrt(np.random.rand(N))
mid_phi = 2 * np.pi * np.random.rand(N)
# the cartesian coordinates of midpoint
mid_x = [mid_radius[i] * np.cos(mid_phi[i]) for i in range(N)]
mid_y = [mid_radius[i] * np.sin(mid_phi[i]) for i in range(N)]
mid_y_square = np.square(mid_radius)
lengths = 2 * np.sqrt(R_sq - mid_y_square)
angle_a = mid_phi + np.arccos(mid_radius/R) # angle of first endpoint
angle_b = mid_phi - np.arccos(mid_radius/R) # angle of second endpoit
x_of_a = R*np.cos(angle_a)
y_of_a = R*np.sin(angle_a)
sol_C_coord_a = list(zip(x_of_a, y_of_a))
x_of_b = R*np.cos(angle_b)
y_of_b = R*np.sin(angle_b)
sol_C_coord_b = list(zip(x_of_b, y_of_b))
solution_C_pct = 100*sum(lengths > triangle_edge_length)/N
```

In [12]:

```
plot_length_distr_and_chords(lengths, R, solution_C_pct,
triangle_edge_length, sol_C_coord_a, sol_C_coord_b, N)
```