Problem 144: Investigating Multiple Reflections of a Laser Beam

Project Euler Problem 144 Solution

EulerSolve provides an optimized solution for Project Euler Problem 144, Investigating Multiple Reflections of a Laser Beam, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The mirror is the ellipse $4x^2+y^2=100$, except for a very small opening at the top where $|x|\le 0.01$. The beam starts at $(0,10.1)$, enters through that opening, and its first actual impact on the mirror is $(1.4,-9.6)$. The task is to count how many reflections occur before the beam reaches the opening again and escapes. This is a deterministic billiard problem inside an ellipse. Once two consecutive points on the path are known, the next segment is forced by the law of specular reflection, so the mathematics is about deriving a stable update rule for one bounce. Mathematical Approach Let $F(x,y)=4x^2+y^2-100$. Every reflection point lies on the curve $F(x,y)=0$. The key is to express one bounce entirely in terms of the previous point and the current point. The Right State Variables A single point on the ellipse is not enough, because the next bounce depends on the incoming direction. If the last two points are $P_{k-1}$ and $P_k$, then the incoming direction at $P_k$ is $$v_{\mathrm{in}}=P_k-P_{k-1}.$$ So the true state of the system is the ordered pair $(P_{k-1},P_k)$. The implementation stores exactly that pair and updates it after every reflection. The Normal Vector and the Reflection Law For an implicit curve $F(x,y)=0$, the gradient is normal to the curve....

Detailed mathematical approach

Problem Summary

The white cell of this problem, a mirror that delays a laser beam, is the ellipse $4x^2+y^2=100$ with the section $-0.01\le x\le0.01$ at the top missing, for light to pass in and out. The beam starts at $(0,10.1)$, just outside the cell, and first hits the mirror at $(1.4,-9.6)$. At every hit it obeys the law of reflection: the arriving and the leaving beam make the same angle with the normal line, which is perpendicular to the tangent of slope $m=-4x/y$ at $(x,y)$. We must find how many times the beam hits the inner surface of the cell before it leaves through the hole, the hit at $(1.4,-9.6)$ included. ^[1]

The whole solution in one sentence: we follow the beam one straight piece at a time, turning it at each hit with the vector form of the law of reflection and finding the next hit as the second root of a quadratic whose first root, the current hit, is known, until the beam lands in the hole.

Mathematical Approach

We describe one bounce by two formulas and check them by hand, find a quantity that never changes during the flight and explains the pattern of the beam, and discuss why rounding does not change the count.

1. The mirror and its normal

Let $F(x,y)=4x^2+y^2-100$. The mirror is the curve $F=0$ and the inside of the cell is $F\lt0$; dividing by $100$ gives $x^2/25+y^2/100=1$, an ellipse with the semi-axes $5$ along the $x$-axis and $10$ along the $y$-axis. The first hit lies on it, because $4\cdot1.96+92.16=100$, and the start point lies outside, because $F(0,10.1)=2.01\gt0$.

The gradient $\nabla F=(8x,2y)$ lists the partial derivatives of $F$, its rates of change in the $x$- and $y$-directions; it is perpendicular to the curve $F=0$ and points to where $F$ grows, that is, out of the cell. ^[2] This agrees with the statement: a line of slope $m=-4x/y$ has the direction $(y,-4x)$, and the dot product $u\cdot w=u_xw_x+u_yw_y$, zero exactly for perpendicular nonzero vectors, gives $(y,-4x)\cdot(8x,2y)=8xy-8xy=0$. So $n=(8x,2y)$ is a normal vector at $(x,y)$, never zero on the mirror, since only the centre has $8x=2y=0$.

2. The law of reflection as a formula

On a horizontal mirror, a beam with direction $(3,-4)$ leaves with direction $(3,4)$: the part along the mirror stays and the part across it is reversed. A tilted mirror does the same, so let a beam with direction $v$ hit the mirror at a point with normal $n$ and split $v$ into a part along the normal and a part along the mirror,

$$v=\frac{v\cdot n}{n\cdot n}\,n+\Bigl(v-\frac{v\cdot n}{n\cdot n}\,n\Bigr),$$

the second part being perpendicular to $n$, as its dot product with $n$ is $v\cdot n-v\cdot n=0$. Reflection keeps the part along the mirror and reverses the part along the normal:

$$v'=v-2\,\frac{v\cdot n}{n\cdot n}\,n.$$

This is the reflection of $v$ in the line perpendicular to $n$; for a normal of length $1$ it reads $v'=v-2(v\cdot n)\,n$. ^[3] Since the two parts are perpendicular, $|v'|=|v|$ by Pythagoras, and $v'\cdot n=-v\cdot n$. So the reversed arriving beam $-v$ and the leaving beam $v'$ make the same angle with the inward normal $-n$, as both angles have the cosine $(v\cdot n)/(|v|\,|n|)$, and they lie on opposite sides of it, as their parts along the mirror are opposite: this is the law of the statement. An arriving beam moves outwards, $v\cdot n\gt0$ (at the first hit by Section 5, at all later ones by Section 3), so $v'\cdot n\lt0$ and it leaves into the cell. Replacing $n$ by $kn$ with $k\ne0$ does not change $v'$.

3. The next hit: a quadratic with a known root

After the bounce at $P=(x,y)$ the beam runs along the points $P+t\,d$ with $t\gt0$, where $d=(d_x,d_y)$ is the new direction $v'$. It meets the mirror again where

$$4(x+t\,d_x)^2+(y+t\,d_y)^2-100=0.$$

Expanding gives $At^2+Bt+C=0$ with $A=4d_x^2+d_y^2$, $B=8x\,d_x+2y\,d_y$ and $C=F(x,y)=0$, because $P$ lies on the mirror. So $t=0$, the point $P$ itself, is one root, and since the two roots add up to $-B/A$ by Vieta's formulas, ^[4] the other root needs no square root:

$$t^{*}=-\frac{8x\,d_x+2y\,d_y}{4d_x^2+d_y^2}.$$

Its numerator is $-d\cdot n=v\cdot n\gt0$ by Section 2, so $t^{*}\gt0$ and the next hit $P+t^{*}d$ lies ahead of the beam. Since $B=-At^{*}$, we have $F(P+t\,d)=At(t-t^{*})$ with $A\gt0$. This is negative for $0\lt t\lt t^{*}$, so between two hits the beam stays inside the cell and meets the mirror, or the hole, only at the ends of its straight pieces; the derivative $At^{*}\gt0$ at $t^{*}$ shows that it arrives at the next hit moving outwards.

4. The way in

The same idea confirms the entry through the hole. The line from $(0,10.1)$ through the first hit consists of the points $(1.4s,\ 10.1-19.7s)$, and $F=0$ becomes

$$395.93\,s^2-397.94\,s+2.01=0.$$

One root is $s=1$, the first hit, and the product of the roots is $2.01/395.93$, so the other is $s\approx0.0050767$, giving approximately $(0.0071073,\ 9.9999899)$: the beam crosses the ellipse inside the hole, and its first contact with the mirror is indeed $(1.4,-9.6)$.

5. Worked example: the first bounce

The beam arrives at $P_0=(1.4,-9.6)$ from $(0,10.1)$, so $v=(1.4,-19.7)$ and $n=(11.2,-19.2)$. Then $v\cdot n=15.68+378.24=393.92$, $n\cdot n=125.44+368.64=494.08$ and $2\cdot393.92/494.08=1231/772$, so

$$v'=v-\tfrac{1231}{772}\,n\approx(-16.4591,\ 10.9155).$$

As Section 2 promised, $|v'|^2=|v|^2=390.05$, and both beams make an angle of about $26.19^\circ$ with the normal. With $d=v'$, Section 3 gives $t^{*}=393.92/1202.7527\ldots$, about $0.32752$, and the second hit

$$P_1=P_0+t^{*}v'\approx(-3.9905976194,\ -6.0249914989).$$

All of this is exact rational arithmetic: $v'=(-15883/965,\ 21067/1930)$, and $P_1$ has the coordinates $-89392050657/q$ and $-134963831648/q$ with $q=22400667565$, for which $F(P_1)=0$ holds exactly. The next hits are $P_2\approx(0.3246,\ 9.9789)$, $P_3\approx(0.4165,\ -9.9652)$ and $P_4\approx(-4.4734,\ 4.4669)$, handy test values for any implementation.

6. The same bounce with slopes

If the tangent makes the angle $\beta$ with the $x$-axis and the arriving beam the angle $\theta$, the leaving beam makes the angle $2\beta-\theta$, and the difference and double-angle formulas for the tangent ^[5] turn $m=\tan\beta$ and $s=\tan\theta$ into the new slope

$$s'=\tan(2\beta-\theta)=\frac{2m-s\,(1-m^2)}{1-m^2+2ms}.$$

At $P_0$, $m=7/12$ and $s=-197/14$ give $s'=-21067/31766$, the slope of $v'$. But slopes become infinite at $(\pm5,0)$ and for vertical beams; the vector formula has no such special cases, which is why all three programs use it.

7. When the beam leaves

By Section 3 the beam leaves exactly when the end of one of its straight pieces lies in the hole, that is, on the ellipse with $|x|\le0.01$ and $y\gt0$. We count the landing points before the first one that passes this test, starting with $(1.4,-9.6)$; the landing in the hole is no hit, as there is no mirror there. Both parts of the test are needed: twice the beam lands near the bottom with $|x|\le0.01$, at about $(0.0060731,\ -9.9999926)$ and $(-0.0033331,\ -9.9999978)$, where the mirror is intact, and a test of $|x|$ alone would stop at the first of them.

8. Exact fractions or floating point

The formulas of Sections 2 and 3 use only the four basic operations, so every hit has rational coordinates and could be computed exactly. But the fractions grow fast: in lowest terms, the two coordinates of $P_1,\dots,P_5$ share denominators of $11$, $36$, $74$, $126$ and $193$ digits, those of $P_{20}$ one of $2904$ digits, roughly in proportion to the square of the number of bounces. The programs therefore use double precision, which rounds each operation to the nearest representable number, for numbers of normal size like these with a relative error of at most $2^{-53}\approx1.1\cdot10^{-16}$. ^[6] Section 11 discusses how these errors add up.

9. A quantity that never changes

The foci of the ellipse are $F_1=(0,c)$ and $F_2=(0,-c)$ with $c=\sqrt{100-25}=5\sqrt3\approx8.660$, and the normal at any point $P$ of the ellipse bisects the angle $F_1PF_2$. ^[7] For a line, let $\delta_1$ and $\delta_2$ be its distances from $F_1$ and $F_2$.

Claim. The product $\delta_1\delta_2$ is the same for all straight pieces of the beam.

Proof. Let the beam arrive at $P$ along the line $\ell$ and leave along $\ell'$. By Section 2, $\ell'$ is the mirror image of $\ell$ in the tangent at $P$, which for a line through $P$ is also its mirror image in the normal. This mirror swaps the lines $PF_1$ and $PF_2$, because the normal bisects the angle between them. So if $\varphi_1$ and $\varphi_2$ are the angles between $\ell$ and the lines $PF_1$ and $PF_2$, then $\varphi_2$ and $\varphi_1$ are those between $\ell'$ and these lines. The distance from a point $Q$ to a line through $P$ is $|PQ|$ times the sine of the angle between the line and $PQ$, so $\delta_1=|PF_1|\sin\varphi_1$, $\delta_2=|PF_2|\sin\varphi_2$, and

$$\begin{aligned}\delta_1'\delta_2'&=|PF_1|\sin\varphi_2\cdot|PF_2|\sin\varphi_1\\&=\delta_1\delta_2.\end{aligned}$$

For the first piece, on the line $19.7x+1.4y=14.14$, the product is $p=132349/975125\approx0.135725$, and hence it is $p$ for every piece.

10. The thin ellipse in the middle

Write the line of a piece as $ax+by=e$ with $a^2+b^2=1$. Its distances from $(0,\pm c)$ are $|e-bc|$ and $|e+bc|$, so if both foci lie on one side of it, $\delta_1\delta_2=e^2-b^2c^2$. The tangent of the ellipse $x^2/\sigma^2+y^2/\tau^2=1$ at its point $(x_1,y_1)$ is $x_1x/\sigma^2+y_1y/\tau^2=1$, ^[7] so the line touches this ellipse exactly when the point $(a\sigma^2/e,\ b\tau^2/e)$ lies on it, that is, when $a^2\sigma^2+b^2\tau^2=e^2$. For an ellipse with the foci of the mirror, $\tau^2=\sigma^2+c^2$, and the condition becomes $\sigma^2=e^2-b^2c^2=\delta_1\delta_2$.

The first piece meets the $y$-axis at $y=10.1$, above $F_1$, so it leaves both foci on one side. The mirror of Section 9 maps the angle $F_1PF_2$ onto itself, so a piece passes between the foci exactly when the next one does, and none of them does. Hence all pieces touch the same ellipse

$$\frac{x^2}{p}+\frac{y^2}{p+75}=1,$$

with the semi-axes $\sqrt p\approx0.3684$ and $\sqrt{p+75}\approx8.6681$. The beam never enters this thin ellipse around the segment $F_1F_2$, which stays as an empty eye in a drawing of the flight. It is the envelope, or caustic, of the trajectory; in an ellipse, trajectories that pass outside the foci have a confocal ellipse as envelope, and Birkhoff showed that billiards in an ellipse are integrable. ^[8]

11. Rounding errors and the exit decisions

In a chaotic billiard such as the Bunimovich stadium, nearby beams separate exponentially fast. ^[8] In the ellipse, a slightly different beam touches a slightly different caustic and advances around it by a slightly different amount per bounce, so nearby beams should drift apart only at a steady rate. In $80$-digit arithmetic, moving the start point to $(10^{-9},\ 10.1)$ while keeping the first hit shifts each coordinate of the hits by at most $9\cdot10^{-6}$, and the largest shift per block of $50$ bounces grows by about $1.2\cdot10^{-6}$ per block. This explains why errors stay small but proves nothing, as rounding adds new errors of about $10^{-15}$ at every bounce. The exact computation under Complexity and Verification settles it: no computed hit is off by more than $7\cdot10^{-10}$ in a coordinate, and every exit test has a far wider margin. Where $|x|\le0.01$, the ellipse gives $|y|\ge9.99998$, so the sign of $y$ is never in doubt; the beam leaves at $x\approx-0.0098488$, $1.5\cdot10^{-4}$ inside the edge of the hole; and the top landing closest to the hole outside it has $|x|\approx0.0165$.

How the Code Works

In the C++ program ^[9] points and directions are structs Point of two double values. reflect_direction is Section 2 with the unnormalised normal $n=(8x,2y)$: from dot $=v\cdot n$ and norm2 $=n\cdot n$ it returns $v-2(v\cdot n)\,n/(n\cdot n)$, with no square root. next_intersection is Section 3: it returns the point for $t^{*}=$num / den, with num $=-(8x\,d_x+2y\,d_y)$ and den $=4d_x^2+d_y^2$.

main first calls parse_arguments, which accepts only --skip-checkpoints; for any other argument it returns false and main returns exit code $1$. Unless the checkpoint is skipped, run_checkpoints performs the first bounce and compares each coordinate of $P_1$ with $x=-3.990597619361616$ and $y=-6.024991498863843$ through nearly_equal, with the absolute tolerance $10^{-9}$; if it fails, main returns exit code $2$. Then main prints the result of solve, which starts with previous at $(0,10.1)$, current at $(1.4,-9.6)$ and the counter reflections at $0$. Each pass of its loop is one bounce: it reflects the direction current - previous at current, finds the next landing point, increases the counter, moves both points forward and stops when std::abs(current.x) <= 0.01 && current.y > 0.0, the test of Section 7. Since the counter grows before the test, the hit at $(1.4,-9.6)$ counts and the landing in the hole does not. Using the start point for the first direction is harmless: it lies on the line of Section 4.

Python runs the same loop in solve, called at the end of the file, and Java directly in main; both use the variables px, py, cx, cy, print the counter bounces at the end and have no options or checkpoint. They scale the normal to length $1$ with math.sqrt or Math.sqrt and use the unit-normal form, which gives the same direction up to rounding, and take the second root $t=-B/A$ from the coefficients $A$ and $B$ of Section 3 (a_coef and b_coef, in Java a and b). All three compute in IEEE double precision, and no loop limits the number of bounces, which is safe because this beam does leave.

Complexity and Verification

If the beam hits the mirror $H$ times, each program makes $H$ passes through its loop, with about thirty floating-point operations each, plus a square root in Python and Java. The time is therefore $O(H)$ and the memory $O(1)$, as only two points, a direction and a counter are kept; each run took a few hundredths of a second on the machine used here, mostly start-up time. All three print the same count, and the C++ checkpoint passes. Independent checks for this article, outside the published solutions:

An exact computation with GMP fractions takes about two minutes because of the growing denominators of Section 8. It finds the same count and the same decision at every exit test, and the exact exit point has $x=-0.0098487613702\ldots$. In each coordinate, the hits of the C++ program differ from the exact ones by at most $4\cdot10^{-10}$; those of the Python and Java programs differ from an $80$-digit computation by at most $7\cdot10^{-10}$.
Along the computed flight the product $\delta_1\delta_2$ of Section 9 stays equal to $p$ up to a relative error of $5\cdot10^{-12}$. By Section 10 this is the same as saying that every computed piece touches the caustic; it checks the invariant, not the positions of the hits.
Run backwards from the exit point, the flight has the same number of hits and ends at the entry point of Section 4, as reflection is reversible.

Footnotes and References

The references supply the background facts; the bounce formulas, the conserved product, its caustic and all numbers are derived or computed in this article.

Project Euler 144 — Laser Beam Reflections. The official statement: the white cell 4x² + y² = 100 without its section −0.01 ≤ x ≤ +0.01 at the top, the beam from (0.0, 10.1) that first hits (1.4, −9.6), the law of reflection and the slope m = −4x/y; it asks how many times the beam hits the internal surface before exiting.
Wikipedia — Gradient. The gradient points where the function increases most quickly and is orthogonal to its level sets.
Wikipedia — Specular reflection. The angle of reflection equals the angle of incidence, and with a unit normal n the reflected direction is d − 2(n · d)n.
Wikipedia — Vieta's formulas. The roots r₁, r₂ of ax² + bx + c satisfy r₁ + r₂ = −b/a and r₁r₂ = c/a.
Wikipedia — List of trigonometric identities. It gives tan(α ± β) = (tan α ± tan β)/(1 ∓ tan α tan β) and tan 2θ = 2 tan θ/(1 − tan² θ).
Wikipedia — Double-precision floating-point format. The 53-bit significand gives 15 to 17 significant digits; rounding to nearest has a relative error of at most 2⁻⁵³ (the machine epsilon), and subnormal numbers trade precision for range.
Wikipedia — Ellipse. The distance from the centre to a focus is c = √(a² − b²), the normal at a point P bisects the angle between the lines to the foci, and the tangent of x²/a² + y²/b² = 1 at (x₁, y₁) is x₁x/a² + y₁y/b² = 1.
Wikipedia — Dynamical billiards. Birkhoff showed that the elliptic billiard is integrable; trajectories outside the foci have a confocal ellipse as envelope, those between them a hyperbola. The Bunimovich stadium is chaotic, with exponentially diverging orbits.
C++ — Euler144.cpp. This immutable revision, with reflect_direction, next_intersection, solve and the first-bounce checkpoint, is the implementation record for this article.

Problem 144 source code

C++

#include <cmath>
#include <iostream>
#include <string>

namespace {

struct Options {
    bool run_checkpoints = true;
};

bool nearly_equal(const double a, const double b, const double eps = 1e-9) {
    return std::abs(a - b) <= eps;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

struct Point {
    double x = 0.0;
    double y = 0.0;
};

Point reflect_direction(const Point& incoming, const Point& point) {
    const double nx = 8.0 * point.x;
    const double ny = 2.0 * point.y;

    const double dot = incoming.x * nx + incoming.y * ny;
    const double norm2 = nx * nx + ny * ny;

    return {
        incoming.x - 2.0 * dot * nx / norm2,
        incoming.y - 2.0 * dot * ny / norm2,
    };
}

Point next_intersection(const Point& current, const Point& direction) {
    const double num = -(8.0 * current.x * direction.x + 2.0 * current.y * direction.y);
    const double den = 4.0 * direction.x * direction.x + direction.y * direction.y;
    const double t = num / den;

    return {
        current.x + t * direction.x,
        current.y + t * direction.y,
    };
}

int solve() {
    Point previous{0.0, 10.1};
    Point current{1.4, -9.6};

    int reflections = 0;
    while (true) {
        const Point incoming{current.x - previous.x, current.y - previous.y};
        const Point outgoing = reflect_direction(incoming, current);
        const Point next = next_intersection(current, outgoing);

        ++reflections;
        previous = current;
        current = next;

        if (std::abs(current.x) <= 0.01 && current.y > 0.0) {
            break;
        }
    }

    return reflections;
}

bool run_checkpoints() {
    const Point previous{0.0, 10.1};
    const Point current{1.4, -9.6};
    const Point incoming{current.x - previous.x, current.y - previous.y};
    const Point outgoing = reflect_direction(incoming, current);
    const Point next = next_intersection(current, outgoing);

    if (!nearly_equal(next.x, -3.990597619361616, 1e-9) ||
        !nearly_equal(next.y, -6.024991498863843, 1e-9)) {
        std::cerr << "Checkpoint failed for first reflected intersection" << '\n';
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve() << '\n';
    return 0;
}

Python

# Problem 144: Investigating multiple reflections of a laser beam
# Laser bounces inside ellipse 4x²+y²=100 until it exits through opening |x|<=0.01, y>0.

import math

def solve():
    # Ellipse: 4x²+y²=100
    px, py = 0.0, 10.1  # previous point (start)
    cx, cy = 1.4, -9.6  # current point (first hit)
    bounces = 0
    
    while True:
        # incoming direction
        dx, dy = cx - px, cy - py
        # normal to ellipse at (cx, cy): gradient of 4x²+y²-100 = (8x, 2y)
        nx, ny = 8*cx, 2*cy
        nlen = math.sqrt(nx*nx + ny*ny)
        nx, ny = nx/nlen, ny/nlen
        # reflect: d' = d - 2*(d·n)*n
        dot = dx*nx + dy*ny
        rdx = dx - 2*dot*nx
        rdy = dy - 2*dot*ny
        # find next intersection with 4x²+y²=100
        # parametric: (cx+t*rdx, cy+t*rdy) on ellipse
        # 4(cx+t*rdx)²+(cy+t*rdy)²=100
        a_coef = 4*rdx*rdx + rdy*rdy
        b_coef = 2*(4*cx*rdx + cy*rdy)
        # c_coef = 4*cx*cx + cy*cy - 100 = 0 (already on ellipse)
        t = -b_coef / a_coef
        nx2 = cx + t*rdx
        ny2 = cy + t*rdy
        bounces += 1
        px, py = cx, cy
        cx, cy = nx2, ny2
        if abs(cx) <= 0.01 and cy > 0:
            break
    
    print(bounces)

solve()

Java

public class Euler144 {
    public static void main(String[] args) {
        double px = 0, py = 10.1, cx = 1.4, cy = -9.6;
        int bounces = 0;
        while (true) {
            double dx = cx - px, dy = cy - py;
            double nx = 8 * cx, ny = 2 * cy, nl = Math.sqrt(nx * nx + ny * ny);
            nx /= nl;
            ny /= nl;
            double dot = dx * nx + dy * ny;
            double rdx = dx - 2 * dot * nx, rdy = dy - 2 * dot * ny;
            double a = 4 * rdx * rdx + rdy * rdy, b = 2 * (4 * cx * rdx + cy * rdy);
            double t = -b / a;
            px = cx;
            py = cy;
            cx += t * rdx;
            cy += t * rdy;
            bounces++;
            if (Math.abs(cx) <= 0.01 && cy > 0)
                break;
        }
        System.out.println(bounces);
    }
}