Project Euler Problem 613
Dave is doing his homework on the balcony and, preparing a presentation about Pythagorean triangles, has just cut out a
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
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