Project Euler Problem 613

Dave is doing his homework on the balcony and, preparing a presentation about Pythagorean triangles, has just cut out a

Project Euler Problem 613

Solution

Answer: 0.3916721504

Let the right triangle have vertices

  • $A=(0,0)$,
  • $B=(30,0)$,
  • $C=(0,40)$,

so the longest side is the hypotenuse $BC$ of length $50$.

For a point $P=(x,y)$ inside the triangle, if the ant chooses a uniformly random direction, the probability it exits through side $BC$ equals the angular width of directions that hit $BC$, divided by $2\pi$.

That angular width is exactly the angle subtended by segment $BC$ at $P$:

$$p(P)=\frac{\angle BPC}{2\pi}.$$

Hence the desired probability is the average of this quantity over the whole triangle:

$$\Pr(BC) = \frac{1}{\text{Area}(\triangle)} \iint_{\triangle} \frac{\angle BPC}{2\pi} ,dA.$$

The triangle area is

$$\frac{30\cdot 40}{2}=600.$$

For $P=(x,y)$,

$$\angle BPC = \cos^{-1} !\left( \frac{(B-P)\cdot(C-P)} {|B-P||C-P|} \right).$$

Explicitly:

$$B-P=(30-x,-y), \qquad C-P=(-x,40-y).$$

So we numerically evaluate

$$\frac{1}{600} \int_0^{30} \int_0^{40-\frac43x} \frac{ \cos^{-1} \left( \frac{(30-x)(-x)+(-y)(40-y)} {\sqrt{(30-x)^2+y^2}\sqrt{x^2+(40-y)^2}} \right) }{2\pi} ,dy,dx.$$

Python implementation

import mpmath as mp

# High precision
mp.mp.dps = 50

def subtended_angle(x, y):
    # vectors from P=(x,y) to B and C
    bx, by = 30 - x, -y
    cx, cy = -x, 40 - y

    # dot product
    dot = bx * cx + by * cy

    # magnitudes
    norm_b = mp.sqrt(bx * bx + by * by)
    norm_c = mp.sqrt(cx * cx + cy * cy)

    # angle BPC
    cos_theta = dot / (norm_b * norm_c)
    return mp.acos(cos_theta)

def integrand(y, x):
    # probability from point (x,y)
    return subtended_angle(x, y) / (2 * mp.pi)

# upper boundary of the triangle:
# y = 40 - (4/3)x
def integrate_x(x):
    ymax = 40 - (4 * x / 3)
    return mp.quad(lambda y: integrand(y, x), [0, ymax])

# integrate over triangle and divide by area
area = 600
result = mp.quad(lambda x: integrate_x(x), [0, 30]) / area

print(result)
print(round(result, 10))

Code walkthrough

  1. subtended_angle(x, y)

Computes the angle $\angle BPC$ using the dot-product formula:

$$\cos\theta = \frac{u\cdot v}{|u||v|}.$$ 2. integrand(y, x)

Converts the angle into a probability by dividing by $2\pi$. 3. Triangle bounds

The hypotenuse is

$$y = 40-\frac43x,$$

so for each $x\in[0,30]$,

$$0\le y\le 40-\frac43x.$$ 4. Double integration

We numerically integrate over the triangle and divide by the area $600$ to average over all starting positions.

Numerical result

The computation gives:

$$0.391672150408749\ldots$$

Rounded to 10 decimal places:

Answer: 0.3916721504