Detailed mathematical approach

Problem Summary

Place the right triangle at \(C=(0,0)\), \(A=(40,0)\), and \(B=(0,30)\), so the legs have lengths \(40\) and \(30\) and the hypotenuse is the side \(AB\) of length \(50\). An interior starting point \(P=(x,y)\) is chosen uniformly from the triangle, and an initial direction is chosen uniformly from \([0,2\pi)\). The ray travels straight until it first meets the boundary.

The task is to find the probability that the exit side is the hypotenuse \(AB\), rounded to \(10\) decimal places. Because the requested precision is much tighter than a simple Monte Carlo simulation would reliably provide, the solution turns the probability into a deterministic area integral and evaluates that integral numerically with high-order quadrature.

Mathematical Approach

The geometric core is simple: for a fixed interior point, the directions that reach the hypotenuse form exactly one angular interval.

Step 1: Describe the Triangle in Coordinates

The triangle is

$T=\left\{(x,y): x\ge 0,\ y\ge 0,\ \frac{x}{40}+\frac{y}{30}\le 1\right\}.$

Its area is

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

For any interior point \(P=(x,y)\), the two endpoints of the target side are \(A=(40,0)\) and \(B=(0,30)\). Everything reduces to understanding the angle subtended by the segment \(AB\) as seen from \(P\).

Step 2: Conditional Probability at a Fixed Point

Fix \(P\in T\). The directions that eventually hit the side \(AB\) are precisely the rays that lie between the rays \(PA\) and \(PB\). Their total angular measure is

$\theta(P)=\angle APB.$

Because the initial direction is uniform on the full circle, the conditional probability is

$\mathbb{P}(AB\mid P)=\frac{\theta(P)}{2\pi}.$

Directions that hit a vertex have measure zero, so they do not affect the average.

Step 3: Compute the Viewing Angle Robustly

Write the vectors from \(P\) to the endpoints of the hypotenuse as

$a=A-P=(40-x,-y),\qquad b=B-P=(-x,30-y).$

The angle between two vectors can be recovered from their oriented area and dot product via

$\theta(P)=\operatorname{atan2}\!\left(\left\lvert a_x b_y-a_y b_x\right\rvert,\ a_x b_x+a_y b_y\right).$

This is exactly the form used by the implementations because it stays numerically stable when the angle is very small or very close to \(\pi\). In this triangle, the numerator also has a geometric meaning:

$\left\lvert a_x b_y-a_y b_x\right\rvert=1200-30x-40y,$

which is twice the area of triangle \(APB\) and is positive for every interior point.

Step 4: Average over All Starting Points

Now average the conditional probability over the uniformly distributed point \(P\):

$p=\frac{1}{\Delta}\iint_T \frac{\theta(P)}{2\pi}\,dA=\frac{1}{2\pi\Delta}\iint_T \theta(x,y)\,dA.$

Since \(\Delta=600\), this becomes

$p=\frac{1}{1200\pi}\iint_T \theta(x,y)\,dA.$

So the problem is reduced to evaluating a two-dimensional integral over the triangle.

Step 5: Map the Triangle to the Unit Square

The numerical integration becomes much cleaner after the substitution

$x=40u,\qquad y=30(1-u)v,\qquad (u,v)\in[0,1]^2.$

This parametrization sweeps the triangle with vertical segments: when \(u=0\) we are on the left side, and when \(u=1\) we collapse to the vertex \(A\). The Jacobian determinant is

$\left\lvert\frac{\partial(x,y)}{\partial(u,v)}\right\rvert=1200(1-u).$

Therefore

$p=\frac{1}{1200\pi}\int_0^1\int_0^1 \theta\!\left(40u,\,30(1-u)v\right)\,1200(1-u)\,dv\,du.$

The denominator \(1200\pi\) that appears in the code comes directly from this formula.

Step 6: Worked Example at One Interior Point

Take \(P=(10,10)\). Then

$a=(30,-10),\qquad b=(-10,20).$

Compute

$\left\lvert a_x b_y-a_y b_x\right\rvert=\left\lvert 30\cdot 20-(-10)(-10)\right\rvert=500,$

$a_x b_x+a_y b_y=30(-10)+(-10)\cdot 20=-500.$

Hence

$\theta(P)=\operatorname{atan2}(500,-500)=\frac{3\pi}{4},$

so the conditional probability at this point is

$\mathbb{P}(AB\mid P)=\frac{3\pi/4}{2\pi}=\frac{3}{8}.$

This is only a local value; the final answer is obtained by averaging \(\theta(P)\) over all interior points.

Step 7: Replace the Integral by Gauss-Legendre Quadrature

After transforming the domain to \([0,1]^2\), the integral is approximated with a tensor-product Gauss-Legendre rule. If \((u_i,w_i)\) are one-dimensional nodes and weights on \([0,1]\), then

$p\approx \frac{1}{1200\pi}\sum_{i=1}^n\sum_{j=1}^n w_i w_j\,\theta\!\left(40u_i,\,30(1-u_i)u_j\right)\,1200(1-u_i).$

Using a large order such as \(n=1024\) gives the required accuracy for the printed \(10\)-decimal result.

How the Code Works

The C++, Python, and Java implementations all begin by constructing Gauss-Legendre nodes and weights on \([0,1]\). They first work on the standard interval \([-1,1]\), use symmetry so only half of the roots need to be found explicitly, and refine the Legendre roots with Newton iteration. The Legendre polynomial values and derivatives are obtained from the standard recurrence, and the nodes and weights are then rescaled from \([-1,1]\) to \([0,1]\).

Once the one-dimensional rule is available, the implementation forms the two-dimensional tensor product. For each outer quadrature node it computes the corresponding \(x\)-coordinate and the Jacobian factor \(1200(1-u)\). For each inner node it computes the matching \(y\)-coordinate, forms the two vectors from the sample point to the endpoints of the hypotenuse, evaluates the angle with the stable \(\operatorname{atan2}\) formula, and adds the weighted contribution to the running integral.

After the double sum is complete, the implementation divides by \(1200\pi\) and formats the probability to \(10\) decimal places. The C++ implementation additionally compares the \(512\)-point and \(1024\)-point quadrature results as a numerical sanity check, while the Python and Java implementations directly evaluate the \(1024\times 1024\) rule.

Complexity Analysis

For quadrature order \(n\), generating the one-dimensional Gauss-Legendre rule costs \(O(n^2)\) arithmetic overall because each root refinement evaluates degree-\(n\) Legendre recurrences. The dominant step is the tensor-product integration, which performs \(n^2\) angle evaluations, so the total running time is \(O(n^2)\). Memory usage is \(O(n)\), since only the one-dimensional nodes and weights plus a few scalars are stored.

Footnotes and References

1. Problem page: Project Euler 613
2. Gauss-Legendre quadrature: Wikipedia - Gauss-Legendre quadrature
3. Geometric probability: Wikipedia - Geometric probability
4. Jacobian matrix and determinant: Wikipedia - Jacobian matrix and determinant
5. Dot product and angle formulas: Wikipedia - Dot product